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" <[EMAIL PROTECTED]> > 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: [EMAIL PROTECTED] >> [mailto:[EMAIL PROTECTED] 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 <[EMAIL PROTECTED]> 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 >> <[EMAIL PROTECTED]> >>> 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 <[EMAIL PROTECTED]> >>> 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 >>> <[EMAIL PROTECTED]> >>>>> 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 <[EMAIL PROTECTED]> >>>>> 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@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, [EMAIL PROTECTED] 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 ______________________________________________ R-help@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.
