On Jun 19, 2012, at 11:39 , anna freni sterrantino wrote: > > > Hi Duncan and Rlist, > I've notice a different behaviour in the invbinomial > you suggest me and invbinomial in stata. > > invbinomial(50,50, 0.4) > Error in uniroot(function(x) pbinom(k, n, x) - p, c(0, 1)) : > f() values at end points not of opposite sign >> invbinomial(50,50, 0.6) > Error in uniroot(function(x) pbinom(k, n, x) - p, c(0, 1)) : > f() values at end points not of opposite sign > > while stata > gen p3=invbinomial(50,50, 0.4) > > . display p3 > 0 > > > > > . gen p4=invbinomial(50,50, 0.6) > > . display p4 > 0 >
Notice that as per your own statement, k==n is outside the domain of definition in Stata, so it is only out of kindness that it is returning any value at all ;-). The root(!) cause is of course that pbinom(n,n,whatever) is identically 1, so solving for any other value is pointless. Modifying Duncan's code to return 0 if k >= n is left as an exercise for the reader... -Peter D. > Thanks > Cheers > > Anna > > > ________________________________ > Da: Duncan Murdoch <murdoch.dun...@gmail.com> > > Cc: Rcran help <r-help@r-project.org> > Inviato: Giovedì 31 Maggio 2012 15:32 > Oggetto: Re: [R] inverse binomial in R > > On 12-05-31 9:10 AM, anna freni sterrantino wrote: >> Hello! >> I'm having some trouble >> trying to replicate in R a Stata function >> >> invbinomial(n,k,p) >> Domain n: 1 to 1e+17 >> Domain k: 0 to n - 1 >> Domain p: 0 to 1 (exclusive) >> Range: 0 to 1 >> Description: returns the inverse of the cumulative binomial; i.e., >> it >> returns the probability of success on one trial >> such >> that the probability of observing floor(k) or >> fewer >> successes in floor(n) trials is p. >> >> I've found some hints on the web like >> http://rwiki.sciviews.org/doku.php?id=guides:tutorials:regression:table >> >> I tried to replicate using qbinom >> the results obtained in >> >>> invbinomial(10,5, 0.5) >>> .54830584 >> >> but with no success. > > I don't think base R has a function like that, though some contributed > package probably does. If you're writing it yourself you'd need to use > uniroot or some other solver, e.g > > invbinomial <- function(n, k, p) { > uniroot(function(x) pbinom(5, 10, x) - p, c(0, 1)) > } > [[alternative HTML version deleted]] > > ______________________________________________ > 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. -- Peter Dalgaard, Professor Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd....@cbs.dk Priv: pda...@gmail.com ______________________________________________ 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.
