From: Katharine Mullen <kate_at_few.vu.nl>

Date: Sat, 19 Apr 2008 18:00:56 +0200 (CEST)

R-help_at_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. Received on Sat 19 Apr 2008 - 16:05:10 GMT

Date: Sat, 19 Apr 2008 18:00:56 +0200 (CEST)

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 <t.yee_at_auckland.ac.nz> 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)
**> >>>
**> >>>
**> >>VGLM linear loop 1 : loglikelihood = -964.1915
**> >>VGLM linear loop 2 : loglikelihood = -964.1392
**> >>VGLM linear 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_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
**> >>
**> >>attached base packages:
**> >>[1] stats4 splines stats graphics grDevices utils datasets methods
**> >>[9] base
**> >>
**> >>other attached packages:
**> >>[1] VGAM_0.7-6 randomForest_4.5-23
**> >>
**> >>______________________________________________
**> >>R-help_at_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.
**> >>
**> >>
**> >>
**>
**>
*

R-help_at_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. Received on Sat 19 Apr 2008 - 16:05:10 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Sat 19 Apr 2008 - 16:30:30 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*