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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.