# Update samples of a Beta with Bernoulli likelihood to the Beta posterior

So this may seem like a strange question, but I have a good reason for it. Nonetheless, I’ll risk the XY problem and describe what I want to do without explaining why I want to do it.

We know that the Beta distribution is the conjugate prior for the Bernoulli distribution, i.e., if $π \sim Beta(\alpha, \beta)$ and $x|π \sim Bernoulli(π)$, then the posterior \begin{align} P(π | x) & \propto p(x | π) \cdot p(π) \\ & = Bernoulli(x; π) \cdot Beta(π; \alpha, \beta) \\ & = Beta(\alpha + x, \beta + 1 - x). \end{align}

Ok. Now, suppose I have a million samples from $Beta(\alpha,\beta)$ and I get a Bernoulli result $x_1 = 1$.

How would I perturb my million samples from the prior to represent the posterior? I know that posterior is $Beta(\alpha+1, \beta)$, but rather than drawing new samples from that posterior, I’d prefer to just adjust my million samples from the prior with $x_1$ to be a Monte Carlo representation of the posterior.

The algebra of the likelihood update above doesn’t immediately suggest to me what to do with Monte Carlo samples to effectively carry it out.

Further detail In reality, I need to do something a little more involved: here’s my processing chain:

1. Draw a million samples from my prior $Beta(\alpha, \beta)$. Let these samples be $p_i$ for $i=1, \cdots, 10^6$ (or some large number).
2. The complication. Raise each sample to a positive power: $p_i ^ \Delta$, for $\Delta \geq 0$.
3. The Monte Carlo ensemble $\lbrace p_i^\Delta \rbrace_i$ still lie on the interval $(0, 1)$ but are no longer Beta-distributed, so fit it to a new $Beta(\alpha', \beta')$ distribution to these new samples. This is my “new” prior.
4. Obtain a Bernoulli realization $x_1$ and update the “new” prior to get a posterior $Beta(\alpha' + x_1, \beta' + 1 - x_1)$.
5. Obtain a million samples from this posterior $Beta(\alpha' + x_1, \beta' + 1 - x_1)$ and go to step 1, awaiting the next Bernoulli experiment $x_2$.

I’d like to condense steps 3–5 by just applying the Bernoulli likelihood on samples of the prior, and just keeping my Monte Carlo ensemble representing the posterior, and then going to step 2 (\Delta can vary). Then, I won’t ever need to fit my ensemble to a Beta.

Having spelled it out like this though, I now realize that even if you gave me a nice way to update samples of a Beta prior with a Bernoulli likelihood, it might not help my situation because my $p_i^\Delta$ won’t be Beta-distributed. I need some way to propagate an arbitrary Monte Carlo ensemble through a Bernoulli likelihood.

• In your example, it's so easy and fast to draw samples that practically any proposal will work. In practice, many methods of performing such simulations do not offer as much information as knowing the full distribution: for instance, you might only know relative densities without knowing the normalization coefficient. For these reasons I suspect you would get more appropriate answers if you were to describe the actual problem you have so that the answers respect the constraints that actually apply to your problem.
– whuber
Apr 11 '17 at 23:47
• If this is a simulation, just reweight your samples using weights that correspond to the ratio of the prior distribution (from which they were drawn) and the posterior distributions (from which you wish they had been drawn). Apr 12 '17 at 0:13
• @whuber you’re right—I updated the question with complete info and while doing so, saw that what I asked for might not apply perfectly to my situation. Please see the update and let me know if this can still be done. Apr 12 '17 at 0:21
• @NeilG thanks for the note, that sounds promising! I am now reading up on reweighting Monte Carlo ensembles. Apr 12 '17 at 0:25

Let $f(x;\alpha,\beta)$ denote the beta distribution parameterized by $\alpha$ and $\beta$. The aim is to express $f(x;\alpha,\beta)$ in terms of $f(x;\alpha-1,\beta)$ and $f(x;\alpha,\beta-1)$.

$$f(x;\alpha, \beta) = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} x^{\alpha - 1}x^{\beta - 1}$$

It turns out the Gamma function satisfies the property $\Gamma(x+1) = x\Gamma(x)$, so this is easy:

$$f(x;\alpha,\beta) = \frac{(\alpha -1 + \beta)\Gamma(\alpha -1 + \beta)}{(\alpha-1)\Gamma(\alpha-1)\Gamma(\beta)}x^{\alpha-2}x^{\beta-1}x^1$$ $$f(x;\alpha,\beta) = \Big( \frac{\alpha-1+\beta}{\alpha-1} x \Big) f(x;\alpha-1,\beta)$$

Similarly,

$$f(x;\alpha,\beta) = \Big( \frac{\alpha-1+\beta}{\beta-1} (1-x) \Big) f(x;\alpha,\beta-1)$$

To adjust your MCMC samples, you can reweight each sample by $\frac{f(x;\alpha,\beta)}{f(x;\alpha-1,\beta)}$ if you observed an alpha result from your bernoulli.

• Here, you’re relating the density of the prior to the posterior—I’d like to update random samples of the prior to the posterior. Apr 12 '17 at 0:06