Hello,
I am using constrOptim to maximize a likelihood function (the values of my
parameter vector must be between zero and one and must sum up to <=1). I am
getting the error 'initial value in 'vmmin' is not finite'. I've tracked it
down to a problem in the function 'R' defined within the constrOptim function.
It is performing a check on my values, and its not passing the check, and
therefore returning NaN. This is not my initial value, but rather during an
internal iteration.
This is the function defined in constrOptim (the following is code from
debugging it)
debugging in: R(theta, theta.old)
debug: {
ui.theta <- ui %*% theta
gi <- ui.theta - ci
if (any(gi < 0))
return(NaN)
gi.old <- ui %*% theta.old - ci
bar <- sum(gi.old * log(gi) - ui.theta)
if (!is.finite(bar))
bar <- -Inf
f(theta, ...) - mu * bar
}
I see that I am getting the gi values that are just machine error below 0
(-8.881784e-16):
Browse[3]> ui.theta
[,1]
[1,] -1
[2,] 1
Browse[3]> ci
[1] -1 0
Browse[3]> gi
[,1]
[1,] -8.881784e-16
[2,] 1.000000e+00
I'm a bit surprised that constrOptim wouldn't be robust to that. I'm wondering
what I can do to get around this, if anything. Is this a sign that my function
that I defined is not behaving correctly in some way (e.g. should I be handling
machine error in the function I defined?). I would note in my example below, I
am only optimizing over a one-dimensional parameter, for which I understand
optim is not optimal, but my function is general and in other situations it
could be a vector of length greater than 1.
I give an example below, as well as my function and session info. The example
is pathological (x, the successes, are equal to m the trials) but I'm running
simulations, and this is being created in my simulations. When I run a large
set of simulations (over many different parameters) I am actually getting about
half the time an error like this, and so I don't know if it is always a
pathological example like this. This is just an example I found where I know it
failed.
Thanks for your help,
Elizabeth Purdom
Example data:
> x <- c(12, 17, 9, 15, 15, 9, 15, 11, 12, 7)
> m <- c(12,17, 9, 15, 15, 9, 15, 11, 12, 7)
> AMat<-matrix(c(1,0,0,1),byrow=T,nrow=2,ncol=2)
> eventTiming(x, m, history = AMat)
Error message:
Error in optim(theta.old, fun, gradient, control = control, method = method, :
initial value in 'vmmin' is not finite
Enter a frame number, or 0 to exit
1: eventTiming(x, m, history = AMat)
2: eventTiming.R#84: constrOptim(theta = init, f = neglik, grad = neggr, ui =
3: eventTiming.R#84: optim(theta.old, fun, gradient, control = control, method
> sessionInfo()
R version 2.15.0 (2012-03-30)
Platform: x86_64-pc-linux-gnu (64-bit)
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=C LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
loaded via a namespace (and not attached):
[1] tools_2.15.0
My function:
eventTiming<-function(x,m,history,
exactAllele=FALSE,normCont=0,verbose=FALSE,coverageCutoff=1,minMutations=10){
A<-history
if(is.null(dim(A))) stop("'history' should be a matrix of size nCopies
x (nEvents +1)")
nCopies<-nrow(A)
K<-ncol(A)-1
possAlleles<-(1:nCopies)/nCopies*(1-normCont)
nMuts<-length(x)
if(length(m)!=nMuts) stop("x and m must be of same length")
if(coverageCutoff<1) stop("coverageCutoff must be at least 1")
if(any(m<coverageCutoff)){
if(verbose) warning("some mutations have no reads, will be
ignored")
nBelowCutoff<-length(which(m<coverageCutoff))
nZero<-length(which(m==0))
x<-x[which(m>=coverageCutoff)]
m<-m[which(m>=coverageCutoff)]
}
else{
nBelowCutoff<-0
nZero<-0
}
if(length(m)<minMutations){
if(length(normCont)>0){
if(length(possAlleles)!=K+1) stop("invalid length for
possible Allele Frequency values")
if(is.null(names(possAlleles)))
names(possAlleles)<-paste("AlleleFreq",1:length(possAlleles),sep="")
return(list(pi=rep(NA,length=K+1),alleleSet=possAlleles,optimOut=NA,nBelowCutoff=nBelowCutoff,nZero=nZero,nMut=length(m)))
}
else
return(list(pi=rep(NA,length=K+1),alleleSet=NA,optimOut=NA,nBelowCutoff=nBelowCutoff,nZero=nZero,nMut=length(m)))
}
###Make 'weight matrix' of P(X=X[i]|X>0, P=p(k,s)) of size nStages x
nCopies for each x --
###an array that is nMuts x nCopies of the prob of X given a certain
allele
probX<-function(p){
dbinom(x=x,size=m,prob=p)
}
wtMat<-sapply(possAlleles,probX) #nMuts x nAlleles
if(is.null(dim(wtMat))){
if(length(x)==1) wtMat<-matrix(wtMat,nrow=1)
else stop("Programming error -- single row wtMat")
}
###Need to convert so in terms of pi[-0]
transMat<-rbind(rep(-1,K),diag(K))
Aid<-A%*%transMat
numGr<-wtMat %*% Aid #N x K matrix
init<-rep(1/(K+1),length=K) #equal distribution; could perhaps choose
better init value
getProbFromTheta<-function(theta){
return(c(1-sum(theta),theta))
}
neglik<-function(theta){
probDist<-as.vector(A%*%getProbFromTheta(theta)) #should be a
vector of probabilities
if(length(probDist)!=K+1) stop("Programming error -- wrong
length")
if(any(probDist< -1e-10)) stop("Programming error -- negative
probabilities")
if(any(probDist> 1+1e-10)) stop("Programming error -- >1
probabilities")
-sum(log(rowSums(sweep(wtMat,2,probDist,FUN="*"))))
}
#numerator of grad because linear in pi -- same regardless of value of
theta
neggr<-function(theta){
#denominator:
probDist<-as.vector(A%*%getProbFromTheta(theta))
if(length(probDist)!=K+1) stop("Programming error -- wrong
length")
if(any(probDist< -1e-10)) stop("Programming error -- negative
probabilities")
if(any(probDist> 1+1e-10)) stop("Programming error -- >1
probabilities")
denom<-rowSums(sweep(wtMat,2,probDist,FUN="*")) #should be a
vector length equal to N
-sum(sweep(numGr,1,denom,"/"))
}
lengthOfTheta<-K
uiL1<-diag(rep(-1,lengthOfTheta),nrow=lengthOfTheta,ncol=lengthOfTheta)
ciL1<-rep(1,lengthOfTheta)
uiGr0<-diag(rep(1,lengthOfTheta),nrow=lengthOfTheta,ncol=lengthOfTheta)
ciGr0<-rep(0,lengthOfTheta)
ui<-rbind(uiL1,uiGr0)
ci<- c(-ciL1,-ciGr0)
#The feasible region is defined by ui %*% theta - ci >= 0.
out<-constrOptim(theta=init,f = neglik, grad=neggr, ui=ui, ci=ci)
return(list(pi=getProbFromTheta(out$par),alleleSet=possAlleles,optimOut=out,nBelowCutoff=nBelowCutoff,nZero=nZero,nMut=length(m)))
}
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
