[R] binomial simulation

[Lucke, Joseph F]

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
> > > > > unit: p(T+,D+) for icu,
> > > > > and the same for reg.
> > > > > I tried:
> > > > >
> > > > > pdeti <- 0
> > > > >
> > > > > pdetr <- 0
> > > > >
> > > > > picu <- pdeti*.6
> > > > >
> > > > > preg <- pdetr*.6
> > > > >
> > > > > dept <- c("icu","reg")
> > > > >
> > > > > icu <- rbinom(200, 1, .6)
> > > > >
> > > > > reg <- rbinom(300, 1, .6)
> > > > >
> > > > > for(i in 1:300) {
> > > > >
> > > > > if(dept=="icu") pdeti==0.95
> > > > >
> > > > > if (dept=="reg") pdetr== 0.80
> > > > >
> > > > > }
> > > > >
> > > > > print(picu)
> > > > >
> > > > > print(preg)
> > > > >
> > > > > and got 50 warnings:
> > > > >
> > > > > the condition has length > 1 and only the
> first
> > > > > element will be used in: if
> > > > > (dept == "icu") pdeti == 0.95
> > > > > the condition has length > 1 and only the
> first
> > > > > element will be used in: if
> > > > > (dept == "reg") pdetr == 0.8
> > > > >
> > > > > I would appreciate any suggestions,
> > > > >
> > > > > thank you,
> > > > >
> > > > > Sigalit.
> > > > >
> > > > >       [[alternative HTML version deleted]]
> > > > >
> > > > >
> ______________________________________________
> > > > > 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<http://www.r-project.org/pos
ting-guide.html>
> > > > > and provide commented, minimal,
> self-contained,
> 
=== 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.