R/penalty.constraint.EXPLN.R

In lisaamrhein/stochprofML: Stochastic Profiling using Maximum Likelihood Estimation

penalty.constraint.EXPLN <-
function(dataset,parameter,smoothingpar=10^5) {
   # this function should only be called if TY>1
   m <- ncol(dataset)
   TY <- length(parameter)/(m+1)
   mu <- parameter[TY:((m+1)*(TY-1))]
   sigma <- parameter[(m+1)*(TY-1)+1]
   lambda <- parameter[(m+1)*(TY-1)+1+(1:m)] 
      
   pen <- 0
   # build a matrix such that the g.th column contains the mu values for gene g
   mu <- matrix(mu,byrow=T,ncol=m)
   # for all genes 1,....m and all types 1,...,TY-1:
   for (g in 1:m) {
      max.datapoint <- max(dataset[,g])
      for (i in 1:(TY-1)) {
         # mode of lognormal i
         mode.ig <- exp(mu[i,g]-sigma^2)
         if (mode.ig<max.datapoint) {
            # sequence of values where to compare the densities
            x <- seq(mode.ig,max.datapoint,(max.datapoint-mode.ig)/500)
            
            # add mode of population i+1 to this sequence if it is on the right of the mode of population i
            # (only applicable if population i+1 is lognormally distributed)
            if (i+1 < TY) {
               # population i+1 has lognormal distribution
               mode.igplus1 <- exp(mu[i+1,g]-sigma^2)
               if (mode.igplus1>mode.ig) {
                  x <- c(x,mode.igplus1)
                  x <- x[order(x)]
               }            
            }
            
            # density of lognormal i at these values
            d.pop1 <- d.sum.of.lognormals(y=x,mu.vector=mu[i,g],sigma.vector=sigma)
            # density of population i+1 at these values
            if (i+1 < TY) {
               # population i+1 has lognormal distribution
               d.pop2 <- d.sum.of.lognormals(y=x,mu.vector=mu[i+1,g],sigma.vector=sigma)
            }            
            else {
               # population i+1 has exponential distribution
               d.pop2 <- d.sum.of.exp.types(y=x,j=1,lambda=lambda[g],logdens=F)            
            }
            # penalize if d.pop2 > d.pop1
            pen <- pen + sum(pmax(0,d.pop2-d.pop1)^2)
            
            # now check whether population i+1 is peaked; in that case, the function d.sum.of.lognormals 
            # most probably smoothed the density too much
            # (only applicable if population i+1 is lognormally distributed)
            if ((i+1 < TY) && (sum(d.pop2)==0)) {
               if (mode.igplus1>mode.ig) {
                  pen <- pen + max(0,dlnorm(mode.igplus1,mu[i+1,g],sigma)-dlnorm(mode.igplus1,mu[i,g],sigma))^2
               }
            }            
         }
      }
   }
   return(smoothingpar*pen)
}