Tell me more ×
Cross Validated is a question and answer site for statisticians, data analysts, data miners and data visualization experts. It's 100% free, no registration required.
up vote 0 down vote favorite

I have an infinite population with unknown mean of successes and failures. I'm drawing 400 times from the population and get 400 successes. Now I want to generate random estimates for the true mean of the population from which I have drawn the 400 successes. Can anyone tell me which probability distribution I have to use, respectively how the function to use in R would look like?

I thought that

rbeta(1,400+1,400-400+1)

might be the right function, but even if I perform this one 10 000 000 times I never get the result 1 (and the same accounts for rbeta(10000000,0+1,400-0+1) and the result 0). So I ask myself why isn't it possible to have a population which consists only of successes or only of failures?

Many thanks in advance!

share|improve this question
(1) Have you tried making some draws from a Beta distribution whose second parameter is close to $0$ and less than $1$? E.g., look at the output of rbeta(10, 400.9, 0.1). (2) Suppose you draw a sample (with replacement) of size $10^{10}$ independently and at random from a population and that all results are successes. (This models the rising of the sun throughout the earth's total potential lifetime.) Are you absolutely certain that the population contains no failures?? – whuber♦ Nov 23 '12 at 16:01
1  
The beta distribution will never generate a value of 1, as it is a continuous distribution and therefore the probability of generating any prespecified value (such as 1) is 0. Are you sure you need to be able to generate an estimate of 1? Or is your real objective to be able to generate a random "sample" of 400 1s based on an estimate derived from the real sample of 400 1s? – jbowman Nov 23 '12 at 16:52
It's unclear what exactly you're trying to do, but my explanation of Bayesian inference on binomial data in this question might be useful. – jerad Nov 24 '12 at 4:30
@jb Please see my comment to Stéphane's answer. The point is that there is an important difference between the Beta distribution as a mathematical object and that distribution as practically applied. In effect, there are Beta distributions that will return zeros in simulations--and that's the setting of this question. – whuber♦ Nov 24 '12 at 23:41
@whuber Ok, I learned, a beta distribution can return zero. – Peter Deplewski Nov 28 '12 at 13:15
show 1 more comment

1 Answer

up vote 1 down vote

I am under the impression your spirit is a little confused. You have to distinguish three things:

  • the real experiment
  • the mathematical (probabilistic/statistical) model for the real experiment
  • the computational aspects of the statistical inference with the mathematical model

For the first point, you claim your population is infinite. So 1) there's no real experiment and your question is purely related to the mathematical model; or 2) you have a real experiment with a huge population and the consideration of an infinite population actually is an assumption of the model (a model can only approximate the real experiment).

You assume a Beta posterior distribution for the proportion of successes. This distribution is a mathematical object, it is a part of the mathematical modeling of your problem. Firstly, a Beta random variable never takes the value $0$ - this is a matter of probability theory; more rigorously we should say that it almost never equals $0$, but the "almost never" probabilistic notion aims to modelize the "never" notion in the real world. Secondly, you use a computer to generate a Beta variable: the numbers you generate are not random, but pseudo-random, and everything is done by the programmer in order that $0$ and $1$ will never be generated (*).

Finally, I suggest you take a look at the probability to sample a success in your infinite population. To do so, you could use the Beta-binomial distribution available in the VGAM package or evaluate this probability using simulations:

> n <- 1000000
> sims <- rep(NA,n)
> for(i in 1:n){
+ p <- rbeta(1,401,1)
+ sims[i] <- rbinom(1,1,p)
+ }
> mean(x)
[1] 0.99754

It is not $1$, why ? Assume you sample $2$ individuals and you get $2$ successes in a population made of $10^6$ individuals. Would you infer that success occurs for every individual of the population ? So what about $400$ successes in an infinite population ?

(*) This point is wrong. See comments below. But for instance a computer never generates the values $0$ and $1$ for a uniform distribution on $[0,1]$.

share|improve this answer
1  
Ask any computer to draw a random double precision value from a $B(10^{-6}, 1-10^{-6})$ distribution: I will give you thousand-to-one odds that it's zero. (It would be a bad bet on your part. :-) – whuber♦ Nov 24 '12 at 23:39
@whuber You're right, except that $10^{-6}$ is not small enough: R returns $0$ when generating $B(10^{-16}, 1-10^{-16})$. – Stéphane Laurent Nov 25 '12 at 7:31
1  
An exact calculation will show that the odds of a double precision underflow (resulting in returning a zero) with $10^{-6}$ are about $1447$ to $1$ in favor. The difference between exponents of $-6$ and $-16$ is important, because in double precision $10^{-16}$ is essentially zero when compared to $1$, so you might be tempted to argue that the problem lies in imprecision in specifying the parameters of the distribution. No such argument is tenable when the parameters are $10^{-6}$ and $1-10^{-6}$. – whuber♦ Nov 25 '12 at 16:36
@whuber: rbeta(1,1e-6,1-1e-6) generates 5.562685e-315 – Stéphane Laurent Nov 25 '12 at 17:00
2  
In double precision, that's a denormalized value considered equal to zero :-). – whuber♦ Nov 25 '12 at 20:53
show 8 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.