# Huso, Dalthorp, Dail and Madsen.2015.Ecological APplications

# Xeq0_Rcode.R

# Calculation of M*:
# M* is based on posterior distribution and alpha.
# The posterior distribution is a function of the prior, g, and x.
# Since we are assuming x = 0, we only need posterior distributions for each combination of g and prior.
# For a given prior, the code below calculates an array with a posterior for each value of g.
# There are three possibilities for g: fixed, uncertain (1.732x), highly uncertain (3x)
#  The degree of uncertainty about g is given in terms of the width of the CI in terms of odds ratios:
#     if Upr and Lwr are upper and lower bounds on CI for g, then odds(Upr)/odds(Lwr) is a measure of uncertainty about g. (NOTE: 1.732 = sqrt(3))

require(tensorA) # tensor algebra

#### initial parameters
# Sometimes need g to range all the way to 1, sometimes needs to stop at 0.99
	g.99<-seq(0.05,0.99, by=.001) # this is assumed to be the same vector throughout the following code.
	g.1<-seq(0.05,1, by=.001) # this is assumed to be the same vector throughout the following code.
	alpha<-c(0.01,0.05,0.1,0.2) # not needed for calculating posteriors but needed for calculating M*
	Mstar<-to.tensor(numeric(length(g.99)*length(alpha)),dim=c(length(g.99),length(alpha))) # array space allocation for later calculations
	names(Mstar)<-c('g','alpha')

#### Defining the uncertainty parameters of g (for g variable only; for fixed g, skip this section)
# uncertain:
	bup.m<-sqrt(3) # "moderate" beta uncertainty parameter
# highly uncertain:
	bup.h<-3 # "highly uncertain" beta uncertainty parameter
	fact<-3

#### fit beta distributions to match the CIs:
# function to define the boundaries of the beta CIs for the g parameter
# and the bounds of the corresponding beta distributions

setBbounds<-function(g,bup,fact){ 
# bup = odds ratio of the upper bound of the 95% interval to the mean (g)
# fact = factor beyond bup to scale the upper bound of the scaled beta distribution to
   ints<-array(dim=c(length(g),2)) #bounds of confidence intervals for each g
   ints[,1]<-g/(g+bup*(1-g)); ints[,2]<-g*bup/(g*bup+1-g)
   phi2<-1/(bup*fact)
   Bu<-g/(g+phi2*(1-g)) # upper bounds on interval that will be scaled & shifted to fit to [0, 1]
   Bl<-phi2*g/(1-g+phi2*g) # lower bounds of fitting intervals
   list(CI=ints,B=cbind(Bl,Bu))
	} # end setBbounds

fitBetas<-function(bounds,g){
   wd<-bounds$B[,2]-bounds$B[,1]; Bl<-bounds$B[,1]
   q025<-(bounds$CI[,1]-bounds$B[,1])/wd; q975<-(bounds$CI[,2]-bounds$B[,1])/wd
   mu<-(g-Bl)/wd; sig2<-((q975-q025)/4)^2 # use mean and variance to determine starting points for optimization
   a<-mu^2/sig2*(1-mu)-mu; b<-a*(1/mu-1); ab<-cbind(a,b) #initial values for optimization
   for (ip in 1:length(g)){
      ab[ip,]<-optim(par=ab[ip,],fn=f1,g=mu[ip],pl=(bounds$CI[ip,1]-Bl[ip])/wd[ip],pu=(bounds$CI[ip,2]-Bl[ip])/wd[ip])$par
	   }
   list(alpha=a,beta=b,CI=bounds$CI,B=bounds$B) #return best fitting parameters for each distribution; NOTE: These are associated with specific bounds
	} # end fitBetas

# find beta parameters that match the quantiles <-- for use in the optimization function for empirical fit to beta distributions
f1<-function(ab,g,pl,pu){
   (qbeta(0.025,shape1=ab[1],shape2=ab[2])-pl)^2+(qbeta(0.975,shape1=ab[1],shape2=ab[2])-pu)^2
	} # end f1

# calculates  P(X = 0 | M = m) = p(X = 0 | g)P(g) assuming g is distributed as a beta RV 
# with given parameters (and scaled to fit in interval [Bl, Bu])
bbin0<-function(g,m,Bl,Bu,a,b){
   (1-g)^m*dbeta((g-Bl)/(Bu-Bl),shape1=a,shape2=b)/(Bu-Bl)
	} # end bbin0

	bounds.h<-setBbounds(g.99,bup.h,fact=fact) ## functions are defined in "utility functions"
	bounds.m<-setBbounds(g.99,bup.m,fact=fact) ## calculations, graphs of g distributions at the end of the file
	betaparms.h<-fitBetas(bounds.h,g.99)
	betaparms.m<-fitBetas(bounds.m,g.99)

calcMstar.gvar<-function(g,priorM,bounds,betaparms,alpha){
   pXe0gM<-to.tensor(numeric(length(priorM)*length(g)),c(m=length(priorM),g=length(g)))
   for (im in 1:dim(priorM)[1]){
      for (ig in 1:length(g)){
         Bl=bounds$B[ig,1]; Bu=bounds$B[ig,2]
         pXe0gM[im,ig]<-integrate(bbin0,lower=Bl,upper=Bu,m=im-1,Bl=Bl,Bu=Bu,betaparms$alpha[ig],b=betaparms$beta[ig])$val
         if (ig%%100==0 & im%%10==0){
            flush.console()
         }
      }
   }
   pMgXe0<-add.tensor(add.tensor(pXe0gM,priorM,op='*'),einstein.tensor(pXe0gM,priorM),op='/') # posteriors for each value of g
   pMltemgXe0<-to.tensor(apply(pMgXe0,F='cumsum',M=2)); names(pMltemgXe0)<-c('m','g')
   Mstar<-to.tensor(numeric(length(g)*length(alpha)),c(g=length(g),alpha=length(alpha)))
   for (ia in 1:length(alpha)){
      Mstar[,ia]<-margin.tensor(pMltemgXe0<1-alpha[ia],'m')
   }
   Mstar
	} # end calcMstar.gvar

##   highly informative : prior distribution on M : initial parameters
#    perfect information about distribution of M~Poisson(3)
	lam<-3
	priorM<-dpois(0:15,lambda=lam)/ppois(15,lambda=lam); priorM<-to.tensor(priorM); names(priorM)<-"m"
	priorM.hi<-priorM
	pXe0gM<-to.tensor(outer(1-g.99,1:dim(priorM)[1]-1,'^')); names(pXe0gM)<-c('g','m')
	pMgXe0<-add.tensor(add.tensor(pXe0gM,priorM,op='*'),einstein.tensor(pXe0gM,priorM),op='/') # posteriors for each value of g
	pMltemgXe0<-to.tensor(apply(pMgXe0,F='cumsum',M=1)); names(pMltemgXe0)<-c('m','g')

	for (ia in 1:length(alpha)){
	   Mstar[,ia]<-margin.tensor(pMltemgXe0<1-alpha[ia],'m')
		}
	Mstar.hi<-Mstar  # M* as a function of g and of alpha

##  moderately informative prior
	priorM<-dnbinom(0:75,size=1.25,p=.1)/pnbinom(75,size=1.25,p=.1); priorM<-to.tensor(priorM); names(priorM)<-"m"
	priorM.mi<-priorM
	pXe0gM<-to.tensor(outer(1-g.99,1:dim(priorM)[1]-1,'^')); names(pXe0gM)<-c('g','m')
	pMgXe0<-add.tensor(add.tensor(pXe0gM,priorM,op='*'),einstein.tensor(pXe0gM,priorM),op='/') # posteriors for each value of g
	pMltemgXe0<-to.tensor(apply(pMgXe0,F='cumsum',M=1)); names(pMltemgXe0)<-c('m','g')
	for (ia in 1:length(alpha)){
	   Mstar[,ia]<-margin.tensor(pMltemgXe0<1-alpha[ia],'m')
		}
	Mstar.mi<-Mstar  # M* as a function of g and of alpha

## variance vs mean
	meanM<-seq(1/20,10,length=200)
	m<-1:30
	pMmat<-array(dim=c(length(meanM),length(m)))

# Poisson
	pMmat.P<-pMmat
	for (i in 1:length(meanM)){
	   pMmat.P[i,]<-dpois(m,meanM[i])
		}
	pMmat<-pMmat.P

# NB 5*mu
	sig2<-5*meanM
	pMmat.NB1<-pMmat
	for (i in 1:length(meanM)){
	   pMmat.NB1[i,]<-dnbinom(m,size=meanM[i]^2/(sig2[i]-meanM[i]),prob=meanM[i]/sig2[i])
		}
	pMmat<-pMmat.NB1

# NB mu + mu^2
	pMmat.NB2<-pMmat
	sig2<-meanM+meanM^2
	for (i in 1:length(meanM)){
	   pMmat.NB2[i,]<-dnbinom(m,size=meanM[i]^2/(sig2[i]-meanM[i]),prob=meanM[i]/sig2[i])
		}
	pMmat<-pMmat.NB2

	priorM<-rep(1/201,201); priorM<-to.tensor(priorM); names(priorM)<-"m"
	priorM.u<-priorM

## fixed g calculations: Mstar
	pXe0gM<-to.tensor(outer(1-g.99,1:dim(priorM)[1]-1,'^')); names(pXe0gM)<-c('g','m')
	pMgXe0<-add.tensor(add.tensor(pXe0gM,priorM,op='*'),einstein.tensor(pXe0gM,priorM),op='/') # posteriors for each value of g
	pMltemgXe0<-to.tensor(apply(pMgXe0,F='cumsum',M=1)); names(pMltemgXe0)<-c('m','g')
	for (ia in 1:length(alpha)){
	   Mstar[,ia]<-margin.tensor(pMltemgXe0<1-alpha[ia],'m')
		}
	Mstar.unif<-Mstar  # M* as a function of g and of alpha

## variable g calculations: Mstar
# moderate uncertainty
	bounds<-bounds.m
	betaparms<-betaparms.m
	Mstar.unif.gm<-calcMstar.gvar(g.99,priorM,bounds,betaparms,alpha)

# high degree of uncertainty
	bounds<-bounds.h
	betaparms<-betaparms.h
	Mstar.unif.gh<-calcMstar.gvar(g.99,priorM,bounds,betaparms,alpha)

## Step 1: Calculate priors for each g (assuming x = 0): result is priorM.pois = an array of prior distributions, one for each g
	M<-0:200
	priorM<-to.tensor(array(numeric(length(g.1)*length(M)),dim=c(length(g.1),length(M)))); names(priorM)<-c('g','m') #one prior for each value of g
	Lx<-seq(.1,100, by=.1) # x-values associated with priorL
	priorL<-1/length(Lx)+numeric(length(Lx)) # what we assume lambda to be before we do any searches
	pXgMpM<-array(dim=c(length(priorL),length(M)))
	for (ig in 1:length(g.1)){
	   for (i in 1:length(M)){
	      pXgMpM[,i]<-dbinom(0,size=i-1,prob=g.1[ig])*dpois(i-1,Lx)
		   }
	   pXgL<-pXgMpM%*%(1+numeric(length(M)))
# this is the posterior of lambda after the first year's search.
#  It looks much like a gamma distribution, which would make a prior M for the next year close to a negative binomial
	   pLgX<-pXgL/sum(pXgL)
	   pMgLpL<-array(dim=c(length(M),length(priorL)))
	   for (i in 1:length(M)){
	      pMgLpL[i,]<-sum(dpois(i-1,Lx)*pLgX)
		   }
	   pM<-pMgLpL%*%(1+numeric(length(priorL)))/length(priorL)
	   priorM[ig,]<-pM/sum(pM)
		}
	priorM.pois<-priorM

## Step 2: calculations of M*: result is Mstar.pois = an array of M* values, one M* for each value of g (and alpha value)
	pXe0gM<-to.tensor(outer(1-g.1,M,'^')); names(pXe0gM)<-c('g','m')
	pMgXe0<-add.tensor(add.tensor(pXe0gM,priorM,op='*'),einstein.tensor(pXe0gM,priorM,only='m',by='g'),op='/')
	pMltemgXe0<-to.tensor(numeric(prod(dim(pMgXe0))),dim=dim(pMgXe0))
	for (ig in 1:length(g.1)){
	   pMltemgXe0[ig,]<-cumsum(pMgXe0[ig,])
		}
	Mstar<-to.tensor(numeric(length(g.1)*length(alpha)),dim=c(length(g.1),length(alpha)))
	names(Mstar)<-c('g','alpha')
	for (ia in 1:length(alpha)){
	   Mstar[,ia]<-margin.tensor(pMltemgXe0<1-alpha[ia],'m')
		}
	Mstar.pois<-Mstar

### same calculation for NB with given 'size' parameter (nb1)
	##   variable p(beta) : overall detection probability : initial parameters
	###   alpha = upper bound of credible interval : initial parameters
	# Calculate priors for each g (assuming x = 0)
	sz<-1 # sz = 1 corresponds to a big variance
	priorM<-to.tensor(array(numeric(length(g.1)*length(M)),dim=c(length(g.1),length(M)))); names(priorM)<-c('g','m') #one prior for each value of g
	Lx<-seq(.1,100, by=.1) # x-values associated with priorL
	priorL<-1/length(Lx)+numeric(length(Lx)) # what we assume lambda to be before we do any searches
	pXgMpM<-array(dim=c(length(priorL),length(M)))
	for (ig in 1:length(g.1)){
	   for (i in 1:length(M)){
	      pXgMpM[,i]<-dbinom(0,size=i-1,prob=g.1[ig])*dnbinom(i-1,mu=Lx,size=sz)
		   }
	   pXgL<-pXgMpM%*%(1+numeric(length(M)))
# this is the posterior of lambda after the first year's search
# it looks just like a gamma distribution, which would make a prior M for the next year a negative binomial
	   pLgX<-pXgL/sum(pXgL)
	   pMgLpL<-array(dim=c(length(M),length(priorL)))
	   for (i in 1:length(M)){
	      pMgLpL[i,]<-sum(dnbinom(i-1,mu=Lx,size=sz)*pLgX)
		   }
	   pM<-pMgLpL%*%(1+numeric(length(priorL)))/length(priorL)
	   priorM[ig,]<-pM/sum(pM)
		}
	priorM.nb1<-priorM #prior distributions assuming M~ negative binomial with mean = lambda and var = mean + mean^2
	pXe0gM<-to.tensor(outer(1-g.1,M,'^')); names(pXe0gM)<-c('g','m')
	pMgXe0<-add.tensor(add.tensor(pXe0gM,priorM,op='*'),einstein.tensor(pXe0gM,priorM,only='m',by='g'),op='/')
	pMltemgXe0<-to.tensor(numeric(prod(dim(pMgXe0))),dim=dim(pMgXe0))
	for (ig in 1:length(g.1)){
	   pMltemgXe0[ig,]<-cumsum(pMgXe0[ig,])
		}
	Mstar<-to.tensor(numeric(length(g.1)*length(alpha)),dim=c(length(g.1),length(alpha)))
	names(Mstar)<-c('g','alpha')
	for (ia in 1:length(alpha)){
	   Mstar[,ia]<-margin.tensor(pMltemgXe0<1-alpha[ia],'m')
		}
	Mstar.nb1<-Mstar

### same calculation for NB with given variance:mean ratio
	# Calculate priors for each g (assuming x = 0)
	mvrat<-5 
	priorM<-to.tensor(array(numeric(length(g.1)*length(M)),dim=c(length(g.1),length(M)))); names(priorM)<-c('g','m') #one prior for each value of g
	Lx<-seq(.1,100, by=.1) # x-values associated with priorL
	priorL<-1/length(Lx)+numeric(length(Lx)) # what we assume lambda to be before we do any searches
	pXgMpM<-array(dim=c(length(priorL),length(M)))
	for (ig in 1:length(g.1)){
	   for (i in 1:length(M)){
	      pXgMpM[,i]<-dbinom(0,size=i-1,prob=g.1[ig])*dnbinom(i-1,size=Lx^2/(mvrat*Lx-Lx),prob=1/mvrat)
		   }
	   pXgL<-pXgMpM%*%(1+numeric(length(M)))
# this is the posterior of lambda after the first year's search
# it looks just like a gamma distribution, which would make a prior M for the next year a negative binomial
	   pLgX<-pXgL/sum(pXgL)
	   pMgLpL<-array(dim=c(length(M),length(priorL)))
	   for (i in 1:length(M)){
	      pMgLpL[i,]<-sum(dnbinom(i-1,size=Lx^2/(mvrat*Lx-Lx),prob=1/mvrat)*pLgX)
		   }
	   pM<-pMgLpL%*%(1+numeric(length(priorL)))/length(priorL)
	   priorM[ig,]<-pM/sum(pM)
		}
	priorM.nb2<-priorM
	pXe0gM<-to.tensor(outer(1-g.1,M,'^')); names(pXe0gM)<-c('g','m')
	pMgXe0<-add.tensor(add.tensor(pXe0gM,priorM,op='*'),einstein.tensor(pXe0gM,priorM,only='m',by='g'),op='/')
	pMltemgXe0<-to.tensor(numeric(prod(dim(pMgXe0))),dim=dim(pMgXe0))
	for (ig in 1:length(g.1)){
	   pMltemgXe0[ig,]<-cumsum(pMgXe0[ig,])
		}
	Mstar<-to.tensor(numeric(length(g.1)*length(alpha)),dim=c(length(g.1),length(alpha)))
	names(Mstar)<-c('g','alpha')
	for (ia in 1:length(alpha)){
	   Mstar[,ia]<-margin.tensor(pMltemgXe0<1-alpha[ia],'m')
		}
	Mstar.nb2<-Mstar