R/bipoilog.R

In poilog: Poisson Lognormal and Bivariate Poisson Lognormal Distribution

    rbipoilog <- function(S,mu1,mu2,sig1,sig2,rho,nu1=1,nu2=1,condS=FALSE,keep0=FALSE){
       
       sim <- function(nr){
         lamx <- rnorm(nr)
         lamy <- rho*lamx+sqrt(1-rho^2)*rnorm(nr)
         x <- rpois(nr,exp(sig1*lamx+mu1+log(nu1)))
         y <- rpois(nr,exp(sig2*lamy+mu2+log(nu2)))
         sel <- x+y 
         xy <- cbind(x,y)
         if (keep0) xy <- cbind(x,y) else xy <- cbind(x[sel>0],y[sel>0])
         return(xy)
       }
      
       
       if (S<1) stop('S is not positive')
       if (!is.finite(S)) stop('S is not finite')
       if ((S/trunc(S))!=1) stop('S is not an integer')
       if (sig1<0) stop('sig1 is not positive')
       if (sig2<0) stop('sig2 is not positive')
       if (nu1<0)    stop('nu1 is not positive')
       if (nu2<0)    stop('nu2 is not positive')
       if (-1>rho | rho>1) stop('rho must be in the range -1 to 1')
       if (!is.logical(keep0)) stop('keep0 should be TRUE or FALSE')
       if (length(c(mu1,mu2,sig1,sig2,rho,nu1,nu2))>7) stop('duplicate values of some parameters')
       
       if (condS){
         simMat <- matrix(NA,0,2)
         fac <- 2
         nr  <- S
         while (nrow(simMat)<S){
           simvals <- sim(nr*fac)
           simMat  <- rbind(simMat,simvals)
           fac     <- (1/(nrow(simvals)/(nr*fac)))*2
           fac     <- ifelse(is.finite(fac),fac,1000)
           nr      <- S-nrow(simvals)
         }
         simMat <- simMat[1:S,]
       }
       
       else simMat <- sim(S)
       return(simMat)
    }
    
       
    rpoilog <- function(S,mu,sig,nu=1,condS=FALSE,keep0=FALSE){
       
       sim <- function(nr){
         lamx <- rnorm(nr)
         x <- rpois(nr,exp(sig*lamx+mu+log(nu)))
         if (!keep0) x <- x[x>0]
         return(x)
       }
       
       if (S<1) stop('S is not positive')
       if (!is.finite(S)) stop('S is not finite')
       if ((S/trunc(S))!=1) stop('S is not an integer')
       if (sig<0) stop('sig is not positive')
       if (nu<0)    stop('nu is not positive')
       
       if (condS){
         simVec <- vector('numeric',0)
         fac <- 2
         nr  <- S
         while (length(simVec)<S){
           simvals <- sim(nr*fac)
           simVec  <- c(simVec,simvals)
           fac     <- (1/(length(simvals)/(nr*fac)))*2
           fac     <- ifelse(is.finite(fac),fac,1000)
           nr      <- S-length(simvals)
         }
         simVec <- simVec[1:S]
       }
       
       else simVec <- sim(S)
       return(simVec)
    }
    
    
   dbipoilog <- function(n1,n2,mu1,mu2,sig1,sig2,rho){
      
      if (length(n1)!=length(n2)) stop('n1 and n2 have unequal length')
      if (any((n1[n1!=0]/trunc(n1[n1!=0]))!=1)) stop('all values of n1 must be integers')
      if (any((n2[n2!=0]/trunc(n2[n2!=0]))!=1)) stop('all values of n2 must be integers')
      if (any(n1<0)) stop('one or several values of n1 are negative')
      if (any(n2<0)) stop('one or several values of n2 are negative')
      
      if (length(c(mu1,mu2,sig1,sig2,rho))>5) stop('vectorization of mu, sig and rho is currently not implemented')
      if (!all(is.finite(c(mu1,mu2,sig1,sig2,rho)))) stop('all parameters should be finite')
      if (-1>rho | rho>1) stop('rho must be in the range -1 to 1')
      if (sig1<0 | sig2<=0) stop('sig1 and/or sig2 is not larger than 0')
      
     .C('poilog2',n1=as.integer(n1),n2=as.integer(n2),mu1=as.double(mu1),mu2=as.double(mu2),
        sig1=as.double(sig1^2),sig2=as.double(sig2^2),rho=as.double(rho),nrN=as.integer(length(n1)),
        val=double(length(n1)), PACKAGE = "poilog")$val
    }
   
   dpoilog <- function(n,mu,sig){
     if (length(mu)>1 | length(sig)>1) stop('vectorization of mu and sig is currently not implemented') 
     if (any((n[n!=0]/trunc(n[n!=0]))!=1)) stop('all n must be integers')
     if (any(n<0)) stop('one or several values of n are negative')
     if (!all(is.finite(c(mu,sig)))) stop('all parameters should be finite')
     if (sig<=0) stop('sig is not larger than 0')
     spos <- which(n<8)
     lpos <- which(n>=8)
     val <- rep(NA, length(n))
     f <- function(x, N) dnorm(x, 0, 1) * dpois(N, exp(x*sig + mu))
     if (length(spos)>0){
       vali <- try(sapply(n[spos], function(n) integrate(f, lower=-Inf, upper=Inf, N=n)$value), silent=TRUE)
       if (inherits(vali, "try-error")){
         vali <- rep(1e-300, length(spos))
         for (i in 1:length(spos)){
           tmp <- try(integrate(f, lower=-Inf, upper=Inf, N=n[i])$value, silent=TRUE)
           if (!inherits(tmp, "try-error")){
             vali[i] <- tmp
           }
         }
       }
       valp <- .C("poilog1", n = as.integer(n[spos]), mu = as.double(mu),
                  sig2 = as.double(sig^2), nrN = as.integer(length(n[spos])),
                  val = double(length(n[spos])), PACKAGE = "poilog")$val
       val[spos] <- apply(cbind(vali,valp),1,max)
     }
     if (length(lpos)>0){
       val[lpos] <- .C("poilog1", n = as.integer(n[lpos]), mu = as.double(mu),
                       sig2 = as.double(sig^2), nrN = as.integer(length(n[lpos])),
                       val = double(length(n[lpos])), PACKAGE = "poilog")$val
     }
     val
   }
    
   poilogMLE <- function(n,startVals=c(mu=1,sig=2),nboot=0,zTrunc=TRUE,
                         method='BFGS',control=list(maxit=1000)){
      
      if (is.matrix(n) | (is.data.frame(n))) stop(paste('n has',ncol(n),'colums, supply a vector or use function bipoilogMLE',sep=' ')) 
      if (length(startVals)!=2) stop('length of startVals is not 2') 
      if (startVals[2]<=0) stop('start value of sig2 is not larger than 0')
      startVals[2] <- log(startVals[2])
      un <- unique(n)
      nr <- rep(NA,length(un))
      for (i in 1:length(un)){ nr[i] <- sum(n%in%un[i]) }
      
      lnL <- function(z,nr){
        if (z[2] < (-372)) z[2] <- -372
        if (z[2] >    354) z[2] <-  354
        b <- 0
        if (zTrunc) b <- log(1-dpoilog(0,z[1],exp(z[2])))
        logL <- -sum((log(dpoilog(un,z[1],exp(z[2])))-b)*nr)
        return(logL)
      }
      
      fit <- optim(startVals,lnL,nr=nr,control=control,method=method)
      
      if (fit$convergence!=0){
        if (fit$convergence==1) stop('the iteration limit has been reached!   try different startVals or increase maxit') 
        if (fit$convergence==10) stop('degeneracy of the Nelder Mead simplex ....')
        else stop(paste('unknown error in optimization', fit$message))
      } 
      
      fit$par <- c(as.numeric(fit$par),1-dpoilog(0,fit$par[1],exp(fit$par[2])))
      est <- list('par'=c('mu'=fit$par[1],'sig'=exp(fit$par[2])),'p'=fit$par[3],'logLval'=-fit$value,'gof'=NULL,boot=NULL)
      
      if (nboot>0){
        cat(paste('estimates: mu: ',round(fit$par[1],3),' sig ',round(exp(fit$par[2]),3),sep=''),'\n')
        cat('********     bootstrapping    ********\n')
        bMat <- matrix(NA,nboot,3)
        colnames(bMat) <- c('mu','sig2','logLval')
        count <- 0
        kat  <- seq(0,nboot,by=100)
        bStartVals <- fit$par
        while (count<nboot){
          bfit <- un <- nr <- NA
          # simulations are conditonal on the number of species in the observed data set :
          sim <- rpoilog(length(n),fit$par[1],exp(fit$par[2]),condS=TRUE,keep0=!zTrunc)
          un <- unique(sim)
          nr <- rep(NA,length(un))
          for (i in 1:length(un)){ nr[i] <- sum(sim%in%un[i]) }
          bfit <- try(optim(bStartVals,lnL,nr=nr,control=control,method=method),silent=TRUE)
          if (inherits(bfit,'try-error')){
            count <- count+1
            bMat[count,] <- c(bfit$par[1],exp(bfit$par[2]),-bfit$value)
            if (count%in%kat) cat('   boot',count,'of',nboot,'\n')
          }
        }
        est$boot <- data.frame(bMat)
        est$gof  <- which(sort(c(est$logLval,bMat[,3]))==est$logLval)/nboot
      }  
      return(est)   
   }
   
  
   
   bipoilogMLE <- function(n1,n2=NULL,startVals=c(mu1=1,mu2=1,sig1=2,sig2=2,rho=0.5),
                           nboot=0,zTrunc=TRUE,file=NULL,
                           method='BFGS',control=list(maxit=1000)){
      
      
      if (is.null(n2)){
        if (!is.matrix(n1) & !is.data.frame(n1)) {stop('n2 is missing and n1 is not a matrix or data frame') }
        if (ncol(n1)!=2) stop('n1 should have 2 columns when n2 is not given')
        n2 <- n1[,2]
        n1 <- n1[,1]
      }
      if (!is.null(n2) & (is.matrix(n1) | is.data.frame(n1))) stop('n1 is a matrix and n2 is specified')
      if (length(n2)!=length(n2)) stop('n1 and n2 must be of equal length')
      
      z <- startVals
      if (length(z)!=5) stop('length of startVals is not 5') 
      if (z[3]<=0) stop('start value of sig1 is not larger than 0')
      if (z[4]<=0) stop('start value of sig2 is not larger than 0')
      if (-1>z[5] | z[5]>1) stop('start value of rho must be in the range -1 to 1')
      
      invRhoTrans <- function(w){
         rho <- (1-exp(-1*w))/(1+exp(-1*w))
         if (rho>0.9999) rho <- 0.9999
         if (rho<(-0.9999)) rho <- -0.9999
         return(rho)
      } 
      rhoTrans    <- function(rho) 1*log((1+rho)/(1-rho))
      
      z[3:4] <- log(z[3:4])
      z[5]   <- rhoTrans(z[5])
 
      xy <- cbind(n1,n2)
      un <- unique(xy)
      nr <- rep(NA,nrow(un))
      for (i in 1:nrow(un)){ nr[i] <- sum(apply(xy,1,function(x) x[1]==un[i,1] & x[2]==un[i,2])) }
      
      lnL <- function(z,nr){
        b <- 0
        if (z[3] < (-372)) z[3] <- -372
        if (z[4] < (-372)) z[4] <- -372
        if (z[3] >   354)  z[3] <-  354
        if (z[4] >   354)  z[4] <-  354
        if (zTrunc) b <- log(1-dbipoilog(0,0,z[1],z[2],exp(z[3]),exp(z[4]),invRhoTrans(z[5])))
        logL <- -sum((log(dbipoilog(un[,1],un[,2],z[1],z[2],exp(z[3]),exp(z[4]),invRhoTrans(z[5])))-b)*nr)
        return(logL)
      }
      
      fit <- optim(z,lnL,nr=nr,control=control,method=method)
      
      if (fit$convergence!=0){
        if (fit$convergence==1) stop('the iteration limit has been reached!   try different startVals or increase maxit') 
        if (fit$convergence==10) stop('degeneracy of the Nelder Mead simplex ....')
        else stop(paste('unknown error in optimization', fit$message))
      }
      
      fit$par <- c(as.numeric(fit$par),1-dpoilog(0,fit$par[1],exp(fit$par[3])),1-dpoilog(0,fit$par[2],exp(fit$par[4])))
      est <- list('par'=c('mu1'=fit$par[1],'mu2'=fit$par[2],'sig1'=exp(fit$par[3]),
             'sig2'=exp(fit$par[4]),'rho'=invRhoTrans(fit$par[5])),'p'=c(fit$par[6],fit$par[7]),'logLval'=-fit$value,'gof'=NULL,boot=NULL)
             
      if (nboot>0){
        cat(paste('estimates: mu1:',round(fit$par[1],3), ' mu2:',round(fit$par[2],3),
                    ' sig1:',round(exp(fit$par[3]),3), ' sig2:',round(exp(fit$par[4]),3),
                    ' rho:',  round(invRhoTrans(fit$par[5]),3), sep=''),'\n')
        cat('********     bootstrapping    ********\n')
        bMat <- matrix(NA,nboot,6)
        colnames(bMat) <- c('mu1','mu2','sig1','sig2','rho','logLval')
        count <- 0
        bStartVals <- fit$par
        while (count<nboot){
          bfit <- un <- nr <- NA
          # simulations are conditonal on the number of species in the observed data set :
          sim <- rbipoilog(length(n1),fit$par[1],fit$par[2],exp(fit$par[3]),exp(fit$par[4]),invRhoTrans(fit$par[5]),condS=TRUE,keep0=!zTrunc)
          un <- unique(sim)
          nr <- rep(NA,nrow(un))
          for (i in 1:nrow(un)){ nr[i] <- sum(apply(sim,1,function(x) x[1]==un[i,1] & x[2]==un[i,2])) }
          bfit <- try(optim(bStartVals,lnL,nr=nr,control=control,method=method),silent=TRUE)
          if (inherits(bfit,'try-error')){
            count <- count+1
            bMat[count,] <- c(bfit$par[1],bfit$par[2],exp(bfit$par[3]),exp(bfit$par[4]),invRhoTrans(bfit$par[5]),-bfit$value)
            cat('   boot',count,'of',nboot,': ',c(bfit$par[1],bfit$par[2],exp(bfit$par[3]),
                                         exp(bfit$par[4]),invRhoTrans(bfit$par[5])),'\n')
          }
          if (!is.null(file)) write.table(bMat[1:count,],file,sep='\t')
        }
        est$boot <- data.frame(bMat)
        est$gof  <- which(sort(c(est$logLval,bMat[,6]))==est$logLval)/nboot
      }  
      
      return(est)   
   }