# prior distribution to the binomial distribution probability distributions urn model

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?

-
(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
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

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]$.

-
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
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
In double precision, that's a denormalized value considered equal to zero :-). –  whuber Nov 25 '12 at 20:53