densities

library(overlapping)

## Loading required package: ggplot2

## Loading required package: testthat

## Warning: package 'testthat' was built under R version 4.0.5

library(survey)

## Loading required package: grid

## Loading required package: Matrix

## Loading required package: survival

## 
## Attaching package: 'survey'

## The following object is masked from 'package:graphics':
## 
##     dotchart

Year1

bioclim data for species 1 and 2 in year 1 and suitability weights (ignore the sort)

data11=sort(rnorm(100000,mean=0.5,sd=0.03))
weights11=rep(1,100000)
weights11=weights11/sum(weights11)

data21=sort(rnorm(100000,mean=0.4,sd=0.1))
weights21=rep(1,100000)
weights21=weights21/sum(weights21)

Year2

data12=sort(rnorm(100000,mean=0.7,sd=0.2))
weights12=1:100000
weights12=weights12/sum(weights12)

data22=sort(rnorm(100000,mean=0.8,sd=0.02))
weights22=1:100000
weights22=weights22/sum(weights22)

data range

lower=min(data11,data12,data22,data21)
upper=max(data11,data12,data22,data21)

Effective population size in each year

ne1=56000
ne2=134000
w1=ne1/(ne1+ne2)
w2=1-w1

Average densities across the two years, weighted by Ne

d11<-density(data11,weights=weights11,n=2048,from=lower,to=upper)
d12<-density(data12,weights=weights12,n=2048,from=lower,to=upper)
d1y<-d11$y*w1+d12$y*w2
d1x<-d11$x
d1<-data.frame(d1x,d1y) 
colnames(d1)=c("x","y")

d21<-density(data21,weights=weights21,n=2048,from=lower,to=upper)
d22<-density(data22,weights=weights22,n=2048,from=lower,to=upper)
d2y<-d21$y*w1+d22$y*w2
d2x<-d21$x
d2<-data.frame(d2x,d2y) 
colnames(d2)=c("x","y")

check average is working

plot(d11,col="red",lwd=3,main="species1")
lines(d12,col="yellow",lwd=3)
lines(d1,col="orange",lwd=3)
legend("topleft",legend=c("year1","year2","avg"),lwd=3,col=c("red","yellow","orange"))

plot(d22,col="blue",lwd=3,main="species2")
lines(d21,col="yellow",lwd=3)
lines(d2,col="green",lwd=3)
legend("topleft",legend=c("year1","year2","avg"),lwd=3,col=c("blue","yellow","green"))

Density functions

dens1<-approxfun(d1)
dens2<-approxfun(d2)

Plot

plot(dens1,xlim=c(lower,upper),ylab="density",col="orange",lwd=2)
curve(dens2,add=T,lwd=2,col="green")
curve(pmin( dens1(x), dens2(x)),add=T,col="black",lwd=2)

Calculate area overlap

integrate(function(x) pmin( dens1(x), dens2(x)),lower,upper )$value/2

## [1] 0.1283371