[R] binomial dist: obtaining probability of success on each trial

David Winsemius dwinsemius at comcast.net
Thu Jan 27 03:51:55 CET 2011


On Jan 26, 2011, at 8:07 PM, Folkes, Michael wrote:

> I'm trying to fathom how to answer two example problems (3.3.2 &  
> 3.3.3) in:
> Krishnamoorthy. 2006. "handbook of statistical distributions with  
> applications"
>
> The first requires calculating single trial probability of success  
> for a binomial distribution when we know:
> trial size=20, successes k=4, P(x<=k)=0.7
> Appreciably all the binomial functions are requiring "prob", which  
> I'm trying to estimate.

Yi might try via successive approximation. Here's a start.

 > dbinom(0:20, 20, prob=0.7)
  [1] 3.486784e-11 1.627166e-09 3.606885e-08 5.049639e-07 5.007558e-06  
3.738977e-05 2.181070e-04 1.017833e-03
  [9] 3.859282e-03 1.200665e-02 3.081708e-02 6.536957e-02 1.143967e-01  
1.642620e-01 1.916390e-01 1.788631e-01
[17] 1.304210e-01 7.160367e-02 2.784587e-02 6.839337e-03 7.979227e-04
 > sum(dbinom(0:20, 20, prob=0.7))
[1] 1
 > sum(dbinom(0:4, 20, prob=0.7))
[1] 5.550253e-06
 > sum(dbinom(0:4, 20, prob=0.2))
[1] 0.6296483
 > sum(dbinom(0:4, 20, prob=0.2))
[1] 0.6296483
 > sum(dbinom(0:4, 20, prob=0.18))
[1] 0.7151181

So you have an interval in which the answer lies and  should be able  
to set up a call to optimize or other solving algorithm.

 > fnbinom <- function(x) abs( sum(dbinom(0:4, 20, prob=x)) - 0.7)
 > optimize(fnbinom, interval=c(0.15, 0.2))
$minimum
[1] 0.1836066

$objective
[1] 6.063185e-05

>
> The second problem is similar:
> prob=0.2, successes k=6, P(x<=k)=0.4, need to calc # trials ('size'  
> in R).

Same sort of approach ought to work.
>
> I'm sure it'll be obvious once somebody explains, but google and R- 
> help archive searches have failed me.  :(
> thanks in advance!
> Michael
> _______________________________________________________
> Michael Folkes
> Salmon Stock Assessment
> Canadian Dept. of Fisheries & Oceans
> Pacific Biological Station
-- 

David Winsemius, MD
West Hartford, CT



More information about the R-help mailing list