Did you notice that the residual function I gave flipped what your
original example problem used as starting values for munb and size? so you
can sometimes afford to give a bad starting value for these param.
in trying to calculate the gradient L-BFGS-B may try evaluate fn with
parameter values outside of the box constraints; to categorically avoid
problems you'd have to know the max. step size for each parameter. it
seems better to use a try statement to catch problems and impose a big
penalty if the try fails, as in
library(VGAM)
phi <- 0
munb <- 72
size <- 0.7
x <- rzinb(200, phi=phi, size=size, munb=munb)
## x should have some non-zero elements
fit <- vglm(x~1, zinegbinomial,trace=TRUE)
fncZiNB <- function(xx, x){
phi <- xx[1]
munb <- xx[2]
size <- xx[3]
yy <- try(dzinb(x=x,phi=phi, size=size, munb=munb,log=TRUE))
if(class(yy) == "try-error" || !all(is.finite(yy))) 1e10 else -sum(yy)
}
# size, the 3rd param must be strictly > 0
solut <- optim(fn=fncZiNB,par=c(0, 72, .7),x=x,method =
"L-BFGS-B", lower=c(0,0,1e-10),upper=c(1,Inf,Inf),
control=list(trace=TRUE),hessian=TRUE)
##optim
solut$par
##vglm
Coef(fit)
On Sat, 19 Apr 2008, Tim Howard wrote:
> Thank you very much for the help and the great suggestion on cleaning up
> the residual function.
>
> Setting the bounds and being very careful with the starting position (on
> my real data) do seem to do the trick, which I wouldn't have guessed
> given the error message that was getting thrown.
>
> Best,
> Tim
>
> >>> Thomas Yee <[EMAIL PROTECTED]> 4/18/2008 6:02 PM >>>
> Hi,
>
> using something like
>
> lower=c(1e-5,1e-5,1e-5),upper=c(1-1e-5,Inf,Inf)
>
> is even more safer since it sometime evaluates the function slightly
> outside the lower and upper limits.
>
> cheers
>
> Thomas
>
>
>
> Katharine Mullen wrote:
>
> >I'm not going to look into what's happening in-depth but it looks like the
> >bounds for your parameters need to be set with care; the below (with
> >slight re-def. of your residual function) gives results that seem to match
> >vglm approximately, with 1e-10 standing in for a bound of 0 and 1-e-10
> >standing in for a bound of 1.
> >
> >library(VGAM)
> >phi <- 0.09
> >size <- 0.7
> >munb <- 72
> >x <- rzinb(200, phi=phi, size=size, munb=munb)
> >
> >fit <- vglm(x~1, zinegbinomial,trace=TRUE)
> >
> >fncZiNB <- function(xx, x){
> > phi <- xx[1]
> > size <- xx[2]
> > munb <- xx[3]
> > -sum(dzinb(x=x,phi=phi, size=size, munb=munb,log=TRUE))
> >}
> >
> > solut <- optim(fn=fncZiNB,par=c(0.09, 72, 0.7),x=x,method =
> >"L-BFGS-B", lower=c(1e-10,1e-10,1e-10),upper=c(1-1e-10,Inf,Inf),
> >control=list(trace=TRUE),hessian=TRUE)
> >
> >##optim
> >solut$par
> >##vglm
> >Coef(fit)
> >
> >On Fri, 18 Apr 2008, Tim Howard wrote:
> >
> >
> >
> >>Dear R-help gurus (and T.Yee, the VGAM maintainer) -
> >>I've been banging my head against the keyboard for too long now, hopefully
> >>someone can pick up on the errors of my ways...
> >>
> >>I am trying to use optim to fit a zero-inflated negative binomial
> >>distribution. No matter what I try I can't get optim to recognize my
> >>initial parameters. I think the problem is that dnbinom allows either
> >>'prob' or 'mu' and I am trying to pass it 'mu'. I've tried a wrapper
> >>function but it still hangs. An example:
> >>
> >>#set up dummy data,load VGAM, which is where I am getting the zero-inflated
> >>neg binom
> >>#I get the same problem without VGAM, using dnbinom in the wrapper instead.
> >>library(VGAM)
> >>phi <- 0.09
> >>size <- 0.7
> >>munb <- 72
> >>x <- rzinb(200, phi=phi, size=size, munb=munb)
> >>
> >>#VGAM can fit using vglm... but I'd like to try optim
> >>
> >>
> >>>fit <- vglm(x~1, zinegbinomial,trace=TRUE)
> >>>
> >>>
> >>VGLMlinear loop 1 : loglikelihood = -964.1915
> >>VGLMlinear loop 2 : loglikelihood = -964.1392
> >>VGLMlinear loop 3 : loglikelihood = -964.1392
> >>
> >>
> >>>Coef(fit)
> >>>
> >>>
> >> phi munb k
> >> 0.1200317 62.8882874 0.8179183
> >>
> >>
> >>
> >>>#build my wrapper function for dzinb
> >>>fncZiNB <- function(phi, size, munb){
> >>>
> >>>
> >>+ -sum(dzinb(x=x,phi=phi, size=size, munb=munb,log=TRUE))}
> >>
> >>
> >>>#test the wrapper, seems to work.
> >>>fncZiNB(phi=0.09, size=0.7, munb=72)
> >>>
> >>>
> >>[1] 969.1435
> >>
> >>#try it in optim
> >>
> >>
> >>>solut <- optim(fn=fncZiNB,par=c(phi=0.09, munb=72, size=0.7),method =
> >>>"L-BFGS-B", hessian=TRUE)
> >>>
> >>>
> >>Error in dnbinom(x = x, size = size, prob = prob, log = log.arg) :
> >> argument "size" is missing, with no default
> >>
> >>I have the same problem with dnbinom.
> >>I'm sure I'm doing something obvious any suggestions?
> >>Session info below. Thanks in advance.
> >>Tim Howard
> >>New York Natural Heritage Program
> >>
> >>
> >>
> >>>sessionInfo()
> >>>
> >>>
> >>R version 2.6.2 (2008-02-08)
> >>i386-pc-mingw32
> >>
> >>locale:
> >>LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
> >>States.1252;LC_M
