inverse binomial in R

[anna freni sterrantino]

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.
Thank you

Cheers

Anna



Anna Freni Sterrantino
Department of Statistics
University of Bologna, Italy
via Belle Arti 41, 40124 BO.
        [[alternative HTML version deleted]]


______________________________________________
[hidden email] 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.

[Duncan Murdoch-2]

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))
}

______________________________________________
[hidden email] 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.

[anna freni sterrantino]

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

Thanks
Cheers 

Anna


________________________________
 Da: Duncan Murdoch <[hidden email]>

Cc: Rcran help <[hidden email]>
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]]


______________________________________________
[hidden email] 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-2]

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 <[hidden email]>
>
> Cc: Rcran help <[hidden email]>
> 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]]
>
> ______________________________________________
> [hidden email] 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: [phone]
Email: [hidden email]  Priv: [hidden email]

______________________________________________
[hidden email] 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.

[Eik Vettorazzi-2]

Hi Anna,
you are at the upper end of the support of the respective binomial
distribution when asking for invbinomial(n,n,p). The probability of
seeing n or less successes in n trials is 1, regardless of the
underlying success rate. It cannot be less, so asking for
invbinomial(50,50, 0.6) makes no sense at all - neither does the result
of stata.

You may insert some exception handling in the function by throwing a
warning or error, when this special case occurs - or even return 0 to
cope the erroneous stata result.

cheers

Am 19.06.2012 11:39, schrieb anna freni sterrantino:

>
>
> 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
>
> Thanks
> Cheers
>
> Anna
>
>
> ________________________________
>  Da: Duncan Murdoch <[hidden email]>
>
> Cc: Rcran help <[hidden email]>
> 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]]
>
>
>
>
> ______________________________________________
> [hidden email] 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.

--
Eik Vettorazzi

Department of Medical Biometry and Epidemiology
University Medical Center Hamburg-Eppendorf

Martinistr. 52
20246 Hamburg

T ++49/40/7410-58243
F ++49/40/7410-57790

--
Pflichtangaben gemäß Gesetz über elektronische Handelsregister und Genossenschaftsregister sowie das Unternehmensregister (EHUG):

Universitätsklinikum Hamburg-Eppendorf; Körperschaft des öffentlichen Rechts; Gerichtsstand: Hamburg

Vorstandsmitglieder: Prof. Dr. Guido Sauter (Vertreter des Vorsitzenden), Dr. Alexander Kirstein, Joachim Prölß, Prof. Dr. Dr. Uwe Koch-Gromus

______________________________________________
[hidden email] 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.

[Nordlund, Dan (DSHS/RDA)]

> -----Original Message-----
> From: [hidden email] [mailto:r-help-bounces@r-
> project.org] On Behalf Of anna freni sterrantino
> Sent: Tuesday, June 19, 2012 2:40 AM
> To: Duncan Murdoch
> Cc: Rcran help
> Subject: Re: [R] inverse binomial in R
>
>
>
> 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
>
> Thanks
> Cheers
>
> Anna
>
>
> ________________________________
>  Da: Duncan Murdoch <[hidden email]>
>
> Cc: Rcran help <[hidden email]>
> 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]]

Anna,

Acccording to the excerpt of Stata documentation that you quote, k must be strictly less than n.  So,
Stata apparently does not return an error (or even a warning) when you supply an invalid input, but the R function does.  How do the functions compare when you supply valid inputs?  What do you want the R function to do when k equals n?


Dan

Daniel J. Nordlund
Washington State Department of Social and Health Services
Planning, Performance, and Accountability
Research and Data Analysis Division
Olympia, WA 98504-5204


______________________________________________
[hidden email] 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.

[anna freni sterrantino]

Hello Peter, Duncan, Dan and RList!
How I suspected... I managed to do the 
"homework" and forced invbinomal to return zero ( as in Stata)
once that k==n.
Thank you for all you replies!

Cheers

Anna

________________________________
 Da: peter dalgaard <[hidden email]>

Cc: Duncan Murdoch <[hidden email]>; Rcran help <[hidden email]>
Inviato: Martedì 19 Giugno 2012 16:00
Oggetto: Re: [R] inverse binomial  in R


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 <[hidden email]>
>
> Cc: Rcran help <[hidden email]>
> 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]]
>
> ______________________________________________
> [hidden email] 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: [phone]
Email: [hidden email]  Priv: [hidden email]
        [[alternative HTML version deleted]]


______________________________________________
[hidden email] 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.