[R] binomial simulation

[Moshe Olshansky]

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/posting-guide.html>
> > > > > and provide commented, minimal,
> self-contained,
> 
=== message truncated ===