[R] binomial simulation

[Prof Brian Ripley]

On Wed, 15 Aug 2007, Moshe Olshansky wrote:

> Thank you - I wasn't aware of this function.
> One can even use lchoose which allows really huge
> arguments (more than 2^1000)!

Using dbinom() for binomial probabilities would be even better, 
and that has a log=TRUE argument to return results on natural log scale.

> dbinom(k,N,p,log=TRUE) + dbinom(m,k,q,log=TRUE)
[1] -92.52584
> log(choose(N,k)*p^k*(1-p)^(N-k)*choose(k,m)*q^m*(1-q)^(k-m))
[1] -92.52584

>
> --- "Lucke, Joseph F" <Joseph.F.Lucke at uth.tmc.edu>
> wrote:
>
>> C is an R function for setting contrasts in a
>> factor.  Hence the funky
>> error message.
>> ?C
>>
>> Use choose() for your C(N,k)
>> ?choose
>>
>> choose(200,2)
>> 19900
>>
>> choose(200,100)
>>  9.054851e+58
>>
>> N=200; k=100; m=50; p=.6; q=.95
>>
> choose(N,k)*p^k*(1-p)^(N-k)*choose(k,m)*q^m*(1-q)^(k-m)
>> 6.554505e-41
>>
>> -----Original Message-----
>> From: r-help-bounces at stat.math.ethz.ch
>> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf
>> Of Moshe Olshansky
>> Sent: Wednesday, August 15, 2007 2:06 AM
>> To: sigalit mangut-leiba; r-help
>> Subject: Re: [R] binomial simulation
>>
>> No wonder that you are getting overflow, since
>> gamma(N+1) = n! and 200! > (200/e)^200 > 10^370.
>> There exists another way to compute C(N,k). Let me
>> know if you need this
>> and I will explain to you how this can be done.
>> But do you really need to compute the individual
>> probabilities? May be
>> you need something else and there is no need to
>> compute the individual
>> probabilities?
>>
>> Regards,
>>
>> Moshe.
>>
>> --- sigalit mangut-leiba <smangut at gmail.com> wrote:
>>
>>> Thank you,
>>> I'm trying to run the joint probabilty:
>>>
>>> C(N,k)*p^k*(1-p)^(N-k)*C(k,m)*q^m*(1-q)^(k-m)
>>>
>>> and get the error: Error in C(N, k) : object not
>> interpretable as a
>>> factor
>>>
>>> so I tried the long way:
>>>
>>> gamma(N+1)/(gamma(k+1)*(gamma(N-k)))
>>>
>>> and the same with k, and got the error:
>>>
>>> 1: value out of range in 'gammafn' in: gamma(N +
>> 1)
>>> 2: value out of range in 'gammafn' in: gamma(N -
>> k) ....
>>>
>>> Do you know why it's not working?
>>>
>>> Thanks again,
>>>
>>> Sigalit.
>>>
>>>
>>>
>>> On 8/14/07, Moshe Olshansky
>> <m_olshansky at yahoo.com>
>>> wrote:
>>>>
>>>> As I understand this,
>>>> P(T+ | D-)=1-P(T+ | D+)=0.05
>>>> is the probability not to detect desease for a
>>> person
>>>> at ICU who has the desease. Correct?
>>>>
>>>> What I asked was whether it is possible to
>>> mistakenly
>>>> detect the desease for a person who does not
>> have
>>> it?
>>>>
>>>> Assuming that this is impossible the formula is
>>> below:
>>>>
>>>> If there are N patients, each has a probability
>> p
>>> to
>>>> have the desease (p=0.6 in your case) and q is
>> the probability to
>>>> detect the desease for a person who
>>> has
>>>> it (q = 0.95 for ICU and q = 0.8 for a regular
>>> unit),
>>>> then
>>>>
>>>> P(k have the desease AND m are detected) = P(k
>> have the desease)*P(m
>>
>>>> are detected / k have
>>> the
>>>> desease) =
>>>> C(N,k)*p^k*(1-p)^(N-k)*C(k,m)*q^m*(1-q)^(k-m)
>>>> where C(a,b) is the Binomial coefficient "a
>> above
>>> b" -
>>>> the number of ways to choose b items out of a
>>> (when
>>>> the order does not matter). You of course must
>>> assume
>>>> that N >= k >= m >= 0 (otherwise the probability
>>> is
>>>> 0).
>>>>
>>>> To generate such pairs (k infected and m
>> detected)
>>> you
>>>> can do the following:
>>>>
>>>> k <- rbinom(N,1,p)
>>>> m <- rbinom(k,1,q)
>>>>
>>>> Regards,
>>>>
>>>> Moshe.
>>>>
>>>> --- sigalit mangut-leiba <smangut at gmail.com>
>>> wrote:
>>>>
>>>>> Hi,
>>>>> The probability of false detection is: P(T+ |
>> D-)=1-P(T+ |
>>>>> D+)=0.05.
>>>>> and I want to find the joint probability
>>>>> P(T+,D+)=P(T+|D+)*P(D+)
>>>>> Thank you for your reply,
>>>>> Sigalit.
>>>>>
>>>>>
>>>>> On 8/13/07, Moshe Olshansky
>>> <m_olshansky at yahoo.com>
>>>>> wrote:
>>>>>>
>>>>>> Hi Sigalit,
>>>>>>
>>>>>> Do you want to find the probability P(T+ = t
>>> AND
>>>>> D+ =
>>>>>> d) for all the combinations of t and d (for
>>> ICU
>>>>> and
>>>>>> Reg.)?
>>>>>> Is the probability of false detection (when
>>> there
>>>>> is
>>>>>> no disease) always 0?
>>>>>>
>>>>>> Regards,
>>>>>>
>>>>>> Moshe.
>>>>>>
>>>>>> --- sigalit mangut-leiba <smangut at gmail.com>
>>>>> wrote:
>>>>>>
>>>>>>> hello,
>>>>>>> I asked about this simulation a few days
>>> ago,
>>>>> but
>>>>>>> still i can't get what i
>>>>>>> need.
>>>>>>> I have 2 units: icu and regular. from icu
>> I
>>> want
>>>>> to
>>>>>>> take 200 observations
>>>>>>> from binomial distribution, when
>> probability
>>> for
>>>>>>> disease is: p=0.6.
>>>>>>> from regular I want to take 300
>> observation
>>> with
>>>>> the
>>>>>>> same probability: p=0.6
>>>>>>> .
>>>>>>> the distribution to detect disease when
>>> disease
>>>>>>> occurred- *for someone from
>>>>>>> icu* - is: p(T+ | D+)=0.95.
>>>>>>> the distribution to detect disease when
>>> disease
>>>>>>> occurred- *for someone from
>>>>>>> reg.unit* - is: p(T+ | D+)=0.8.
>>>>>>> I want to compute the joint distribution
>> for
>>>>> each
>>
> === message truncated ===
>
> ______________________________________________
> R-help at stat.math.ethz.ch 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.
>

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595