# This file includes my own functions for diameter distribution modeling. 
# You are welcome to use them but I take no responsibuility about inpropewr functioning or 
# wrong results. The functions come absolute with no warranty. 

# Reynolds error index
# d includes true diameters, cdf is a function returning the unweighted cdf
# power gives the weights (2 gives BA weight, 0 the unweighted),
# dlim the range of diameters.
# cwidth the class width (make sure dlim is a multiple of dlim).
errind<-function(d,cdf=function(x) pweibull(x,5,20),power=0,dlim=c(-20,80),cwidth=2) {
        limits<-seq(dlim[1],dlim[2],cwidth)
        n<-length(limits)-1
        means<-(limits[-1]+limits[-(n+1)])/2
        dp<-d^power
        sumdp<-sum(dp)
        ftrue<-sapply(1:n,function(i) sum(dp[d>=limits[i]&d<limits[i+1]]))/sumdp
        Fpred<-cdf(limits)
        fpred<-Fpred[-1]-Fpred[-n]
        fpred<-fpred*means^power/sum(fpred*means^power)
        sum(abs(fpred-ftrue))
        }

# 2 parameter Weibull -logL
nLLweibull<-function(x, shape=5,scale=20) {
            -sum(dweibull(x,shape=shape,scale=scale,log=TRUE))
            }    
            
# 2 parameter Weibull density with analytic diffeentiation
mydweibull<-deriv(~(shape/scale)*(d/scale)^(shape-1)*exp(-(d/scale)^shape),
                  c("shape","scale"),function(d,shape,scale){})
  
# 3-parametric Weibull density
dweib3<-function(x,shape=5,scale=20,loc=0) {
        dweibull(x-loc,shape,scale)
        }
        
# 3 parameter weibull -logL
nLLweib3<-function(x, shape=5,scale=20,loc=0) {
            -sum(dweibull(x-loc,shape=shape,scale=scale,log=TRUE))
            }   

# Fits two-parameter weibull distribution to tree diameter data using MLE starting from values 5 and 20
# for shape and scale, respectively
fitw2<-function(d) {
       est<-mle(function(shape=5,scale=20) nLLweibull(d,shape,scale))
       if (class(est)=="try-error") list(par=rep(NA,2),neg2LL=NA,conv=NA) 
       else list(par=coef(est),neg2LL=2*attributes(est)$min,conv=attributes(est)$details$convergence)
       }

# Density of the logit-logistic distribution of Rennolls and Wang.           
# sigma and my unbounded but should be close to 0, say between -5 and 5
# ksi and lambda >0, lambda>ksi
# ksi=minimum, lambda=maximum, sigma controls curtosis and myy skewness.
dll<-function(x,lambda,sigma,ksii,myy) {
     sigma<-exp(sigma)/(1+exp(sigma))
     f<-(lambda-ksii)/sigma*
     1/((x-ksii)*(lambda-x))*
     1/(exp(-myy/sigma)*((x-ksii)/(lambda-x))^(1/sigma)+
        exp(myy/sigma)*((x-ksii)/(lambda-x))^(-1/sigma)+2)
     f[x<=ksii|x>=lambda]<-0
     f
     }

# negative log likelihood     
nLLll<-function(x,lambda,sigma,ksii,myy) {
#      cat(log(dll(x,lambda,sigma,ksii,myy)),"\n")
       value<--sum(log(dll(x,lambda,sigma,ksii,myy)))
       if (value==Inf) value<-10^16
       cat(lambda,sigma,ksii,myy,value,"\n")
       value
       }     
         
# derivatives of neg ll
# theta=c(lambda,sigma,ksii,myy)
gradnLLll<-function(x,theta) {
                l<-theta[1]
                s<-theta[2]
                k<-theta[3]
                m<-theta[4]
                awful<-exp(-m/s)*((x-k)/(l-x))^(1/s)+exp(m/s)*((l-x)/(x-k))^(1/s)+2
                c(
                # with respect to lambda
                sum(-1/(l-k)+1/(l-x)+
                1/awful*
                (-exp(-m/s)*(1/s)*((x-k)/(l-x))^((1-s)/s)* (x-k)*(l-x)^(-2)+
                 +exp(m/s)* (1/s)*((x-k)/(l-x))^((-1-s)/s)*(x-k)*(l-x)^(-2))),
                # with respect to sigma
                sum(1/s+
                1/awful*
                (-(exp(-m)*((x-k)/(l-x)))^(1/s)*(-m+log((x-k)/(l-x)))*s^(-2)+
                 -(exp(m)*((l-x)/(x-k)))^(1/s)*(m+log((l-x)/(x-k)))*s^(-2))),
                #with respect to ksii
                sum(1/(l-k)-1/(x-k)+
                1/awful*
                (-exp(-m/s)*(1/s)*((x-k)/(l-x))^((1-s)/s)/(l-x)+
                 exp(m/s)*(1/s)*((x-k)/(l-x))^((-1-s)/s)/(l-x))),
                # with respect to myy
                sum(1/awful*
                (((x-k)/(l-x))^(1/s)*exp(-m/s)*(-1/s)+
                 ((x-k)/(l-x))^(-1/s)*exp(m/s)*(1/s)))
                )
                }

# cdf of logit-logistic distribution     
pll<-function(x,lambda,sigma,ksii,myy) {
     sigma<-exp(sigma)/(1+exp(sigma))
     F<-1/(1+exp(myy/sigma)*((x-ksii)/(lambda-x))^(-1/sigma))
     F[x<=ksii]<-0
     F[x>=lambda]<-1
     F
     }

# a function for fitting logit-logistic distribution to tree diameter data.
fitll<-function(d,start=NA) {
       ll<-function(lambda,sigma,ksii,myy) nLLll(d,lambda,sigma,ksii,myy)
       if (is.na(start)) start<-list(lambda=max(d)+5,sigma=0.5,ksii=max(0,min(d)-2),myy=0)
       grad<-function(theta) gradnLLll(d,theta)
       est<-try(mle(ll,gr=grad,method="L-BFGS-B",start=start,lower=c(max(d)+0.1,0.1,0,-5),upper=c(45,0.999,max(0,min(d)-0.1),5),control=list(maxit=200)))
       if (class(est)=="try-error") list(par=rep(NA,4),neg2LL=NA,conv=NA) 
       else list(par=coef(est),neg2LL=2*attributes(est)$min,conv=attributes(est)$details$convergence)
       }

# Another function for the same purpose
#  sigma is bounded using logit transformation
fitll<-function(d,start=NA) {
       ll<-function(lambda,sigma,ksii,myy) nLLll(d,lambda,sigma,ksii,myy)
       if (is.na(start)) start<-list(lambda=max(d)+5,sigma=0,ksii=max(0,min(d)-2),myy=0)
       est<-try(mle(ll,method="L-BFGS-B",start=start,lower=c(max(d)+0.1,-Inf,0,-5),upper=c(45,3.5,max(0,min(d)-0.1),5),control=list(maxit=200)))
       if (class(est)=="try-error") list(par=rep(NA,4),neg2LL=NA,conv=NA) 
       else list(par=coef(est),neg2LL=2*attributes(est)$min,conv=attributes(est)$details$convergence)
       }

# quantile function of the logit-logistic distribution       
qll<-function(p,lambda,sigma,ksii,myy) {
         apu<-sapply(1:length(p),function (i) updown(l=ksii,u=lambda,
                     fn=function(x) pll(x,lambda,sigma,ksii,myy)-p[i],crit=10)$par)
         apu
         }

# random number generation from the logit-logistic distribution         
rll<-function(n,lambda,sigma,ksii,myy) {
         qll(runif(n),lambda,sigma,ksii,myy)
         }

# the relatinship between N and G when tree diameter follows the logit-logistic distribution
npergll<-function(lambda,sigma,ksii,myy) {
     1/integrate(function(x) pi*x^2/40000*dll(x,lambda,sigma,ksii,myy),ksii,lambda)$value
     }

# Rennolls&Wang parameterization of Johnsons SB distribution
# lamda=maximum, not range
# desity
dsb<-function(x,ksi,lambda,delta,gamma) {
     value<-rep(NA,length(x))
     value[x>ksi&x<lambda]<-1/(delta*sqrt(2*pi))*(lambda-ksi)/((lambda-x[x>ksi&x<lambda])*(x[x>ksi&x<lambda]-ksi))*
     exp(-1/2*((log((x[x>ksi&x<lambda]-ksi)/(lambda-x[x>ksi&x<lambda]))-gamma)/delta)^2)
     value[x<=ksi|x>=lambda]<-0
     value
     }
     
# c.d.f
psb<-function(x,ksi,lambda,delta,gamma) {
     y<-rep(NA,length(x))
     y[x<=ksi]<-0
     y[x>=lambda]<-1
     y[x>ksi&x<lambda]<-sapply(x[x>ksi&x<lambda],
                        function(y) integrate(function(x) dsb(x,ksi,lambda,delta,gamma),ksi,y)$value)
     y
     }
# negative log-likelihood    
nLLsb<-function(x,ksi,lambda,delta,gamma) {
          value<--sum(log(dsb(x,ksi,lambda,delta,gamma)))
          if (value==Inf) value<-10^16
#          cat(ksi,lambda,delta,gamma,value,"\n")
          value
          }

# gradients of the negative log-likelihood         
gradnLLsb<-function(x,theta) {
               k<-theta[1]
               l<-theta[2]
               d<-theta[3]
               g<-theta[4]
               big<-(log((x-k)/(l-x))-g)/d
               c(
               # with respect to ksi
               sum(1/(l-k)-1/(x-k)-big/d/(x-k)),
               # with respect to lambda
               sum(-1/(l-k)+1/(l-x)-big/(d*(l-x))),
               # with respect to delta
               sum(1/d-big*(log((x-k)/(l-x))-g)*d^(-2)),
               # with respect to gamma
               sum(-big/d)
               )
               }

# a function for fitting Johnsons SB distribution to tree diameter data using MLŽ.               
fitsb<-function(d) {
       d<-d
       ll<-function(ksi,lambda,delta,gamma) nLLsb(d,ksi,lambda,delta,gamma)
       start<-list(ksi=max(0,min(d)-2),lambda=max(d)+5,delta=1,gamma=0)
       grad<-function(theta) gradnLLsb(d,theta)
       est<-try(mle(ll,gr=grad,start=start,method="L-BFGS-B",lower=c(0,max(d)+0.1,0.01,-10),upper=c(max(0,min(d)-0.1),Inf,Inf,10),control=list(maxit=500)))
       if (class(est)=="try-error") list(par=rep(NA,4),neg2LL=NA,conv=NA) 
       else list(par=coef(est),neg2LL=2*attributes(est)$min,conv=attributes(est)$details$convergence)
       }
          

resplot<-function(x,res,limits=NA,ylim=range(res),xlab="x",ylab="res",main="",points=TRUE) {
         if (points) {
            plot(x,res,col="red",ylim=ylim,xlab=xlab,ylab=ylab,main=main)
            mywhiskers(x,res,se=FALSE,add=TRUE,lwd=1,limits=limits)
            } else {
            mywhiskers(x,res,se=FALSE,add=FALSE,lwd=1,limits=limits,xlab=xlab,ylab=ylab)
            }
         mywhiskers(x,res,add=TRUE,lwd=5,limits=limits)         
         abline(0,0,col="gray")
         }


        
# Gussian kernel. Same as density() except that x-values are given as input        
dgkernel<-function(x,d,bw=bw.nrd0(d)) {
          n<-length(d)
          sapply(x,function(a) 1/(n*bw)*sum(dnorm((a-d)/bw)))
          }

# c.d.f. of a gaussian kernel          
pgkernel<-function(x,d,bw=bw.nrd0(d)) {
          n<-length(d)
          sapply(x,function(a) 1/n*sum(pnorm((a-d)/bw)))
          }
          
# quantile function of a gaussian kernel    
qgkernel<-function(p,d,bw=bw.nrd0(d)) {
         apu<-sapply(1:length(p),function(i) updown(l=0,u=60,
                     fn=function(x) pgkernel(x,d,bw)-p[i],crit=10)$par)
         apu
         }