2
$\begingroup$

The Bayesian formula is given as the following simple way.

$${\mathsf P}(a\mid x) ~=~ \dfrac{{\mathsf P}(x\mid a)~{\mathsf P}(a)}{{\mathsf P}(x)}$$

A factory makes pencils.
prior probability: defective pencils manufactured by the factory is 30%.
To check 10 pencils ,2 defective pencil found.

a is event : defective rate of pencils.
x is sample to check the pencils.
prior probability : P(a) = 0.3
P(x|a) follows binomial distribution, expressed in R.

n <- 10
x <- 2
choose(n,x)*0.3^2*(1-0.3)^(10-2)
0.2334744

P(x|a) =0.233

How can get P(x) to calculate the posterior probability ${\mathsf P}(a \mid x)$ in this example with the Bayesian formula?

$\endgroup$
  • 2
    $\begingroup$ $a$ is either an event or a parameter, i.e. the probability of a pencil being defective. When it is the parameter, as in your usage in $P(X=x|a)$, we should be talking about the prior distribution of $a$, i.e. $p_a(a')$. $\endgroup$ – gunes Jul 23 at 6:07
1
$\begingroup$

In order to treat this situation as a problem in Bayesian inference, the probability $\theta = P(\text{Defective})$ must be considered as a random variable. Its prior distribution cannot be taken as degenerate with $P(\theta=0.3)=1.$ In that case, binomial data could not be used to modify the prior distribution, in order to obtain a posterior distribution. (I believe that is the crux of the comment by @gunes.) You could say that the constant $\theta = 0.3$ leads to the probability $$P(X = 2\,|\, n=10,\theta=.3) =0.2334,$$ computed using R as in your Question or more conveniently, because dbinom is a binomial PDF in R, as below. But that is a basic probability computation, not Bayesian inference.

dbinom(2, 10, .3)
[1] 0.2334744

Useful prior distributions. There are several ways to specify a prior distribution for $\theta$ in order to do Bayesian inference, which are consistent with your idea that, roughly $\theta \approx 0.3.$

Beta prior. For example, the distribution $\mathsf{Beta}(3, 7)$ has density function $p(\theta) = K\theta^{3-1}(1-\theta)^{7-1}$, for $0 < \theta < 1 ,$ So $E(\theta) = 3/(3+7) = 0.3$ and the distribution and places 95% of its probability "near" 0.3 [in the interval $(0.075, 0.600)].$ (For information about beta distributions see Wikipedia.)

qbeta(c(.025,.975), 3, 7)
[1] 0.07485463 0.60009357

Binomial likelihood. Now, if you observe $X = 2$ defective pencils out of $n = 10,$ that corresponds to the binomial likelihood $$p(x\,|\,\theta) \propto \theta^2(1-\theta)^8,$$ where the symbol $\propto$ indicates that we have omitted the norming constant ${10 \choose 2},$ thus specifying the likelihood function only 'up to a constant'.

Beta posterior. Then, by Bayes' Theorem, the posterior distribution is proportional to the product of the prior beta density $p(\theta)$ and the binomial likelihood function.

$$p(\theta\, |\, x) \propto p(\theta)\times p(x\, |\, \theta) \propto \theta^{3-1}(1\theta)^{7-1} \times \theta^2 (1-\theta)^8 \\ \propto \theta^{5-1}(1-\theta)^{15-1},$$ where we recognize the final term as the 'kernel' (PDF without norming constant) of the distribution $\mathsf{Beta}(5, 15),$ which has mean $5/(5+15) = 0.25$ and places 95% of its probability in the interval $(0.091,0.456).$ This is sometimes called a 'Bayesian 95% posterior probability interval' for $\theta.$

qbeta(c(.025, .975), 5, 15)
[1] 0.09146578 0.45565308

Taken together, the information provided by the prior distribution and the information provided by the data have given a posterior distribution that concentrates its probability closer to 0.25 than to 0.30.

If you want a prior distribution with probability more concentrated near $\theta = 0.3,$ then you might choose the prior distribution $\mathsf{Beta}(15,35)$ with $E(\theta) = 0.3,$ and 95% of its probability in $(0.183,0.433).$ Then with $X=2$ successes out of $n=10$ trials, the posterior mean would be $E(\theta) = 0.2833$ and the 95% posterior probability interval is $(0.178, 0.403).$ The 'stronger' prior has a greater influence on the posterior distribution.

qbeta(c(.025,.975), 15, 35)
[1] 0.1825194 0.4326297
qbeta(c(.025,.975), 15+2, 35+8)
[1] 0.1776482 0.4026935

Simulating posterior information. It is easy to find the posterior distribution in this situation because the beta and binomial distributions are 'conjugate' (have compatible mathematical forms), so that we do not need to compute the denominator in Bayes' Theorem.

The following simulation in R approximates the joint distribution if $\theta$ and $X,$ and then finds the posterior distribution given that $X = 2.$ With a million iterations, results are accurate to a couple of decimal places.

set.seed(723)     # for reproducibility
m = 10^6; th = x = numeric(m)
for(i in 1:m) {
   th[i] = rbeta(1, 3, 7)
   x[i] = rbinom(1,10,th[i]) }
mean(th[x==2]);  5/20    
[1] 0.250043    # aprx posterior mean = 0.25
[1] 0.25        # exact posterior mean
quantile(th[x==2], c(.025, .975))
      2.5%      97.5% 
0.09101452 0.45530780 # aprx 95% posterior prob int
qbeta(c(.025,.975), 3+2, 7+8)
[1] 0.09146578 0.45565308 # exact 95% PI

A discrete prior distribution. You could pick a prior distribution on $\theta$ that takes values $(.2, .3, .4)$ with respective probabilities $(.25, .50, .25).$ Then the prior mean is $E(\theta) = 0.3.$ Finding the posterior distribution requires some simple, but perhaps tedious, computation using the elementary (discrete) version of Bayes' Theorem. The following simulation approximates the posterior distribution, given $X = 2.$

set.seed(2019)
th.val = c(.2,.3,.4); prior=c(.25,.50,.25)
m = 10^6;  th = x = numeric(m)
for(i in 1:m) {
   th[i] = sample(th.val,1,pr=prior)
   x[i] = rbinom(1, 10, th[i]) }
table(th[x==2])/sum(x==2)

      0.2       0.3       0.4 
0.3409187 0.5235124 0.1355689 
mean(th[x==2])
[1] 0.279465

A weaker prior with probabilities $(1/3,1/3,1/3)$ gives results:

table(th[x==2])/sum(x==2)

      0.2       0.3       0.4 
0.4616765 0.3548434 0.1834800 
mean(th[x==2])
[1] 0.2721803

Once again, what you cannot do, if you want actual Bayesian inference is to make $P(\theta = 0.3) = 1.$ Then the data would be completely irrelevant.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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