# Posterior sampling without using pm.Potential in pyMC3

I'm going through the Price Is Right example in chapter 5 of Probabilistic Programming & Bayesian Methods for Hackers and I have problems understanding the solution. I have tried to change the implementation of the solution in a way that I thought that should be also correct but I got results that deviate from those in the script.

The description of the example is given below:

Example: Optimizing for the Showcase on The Price is Right

Bless you if you are ever chosen as a contestant on the Price is Right, for here we will show you how to optimize your final price on the Showcase. For those who forget the rules:

1. Two contestants compete in The Showcase.
2. Each contestant is shown a unique suite of prizes.
3. After the viewing, the contestants are asked to bid on the price for their unique suite of prizes.
4. If a bid price is over the actual price, the bid's owner is disqualified from winning.
5. If a bid price is under the true price by less than \\$250, the winner is awarded both prizes.

The difficulty in the game is balancing your uncertainty in the prices, keeping your bid low enough so as to not bid over, and trying to bid close to the price.

Suppose we have recorded the Showcases from previous The Price is Right episodes and have prior beliefs about what distribution the true price follows. For simplicity, suppose it follows a Normal:

$$\rm{True Price} \sim \text{Normal}(\mu_p, \sigma_p )$$

In a later chapter, we will actually use real Price is Right Showcase data to form the historical prior, but this requires some advanced PyMC3 use so we will not use it here. For now, we will assume $$\mu_p = 35 000$$ and $$\sigma_p = 7500$$.

We need a model of how we should be playing the Showcase. For each prize in the prize suite, we have an idea of what it might cost, but this guess could differ significantly from the true price. (Couple this with increased pressure being onstage and you can see why some bids are so wildly off). Let's suppose your beliefs about the prices of prizes also follow Normal distributions:

$$\text{Prize}_i \sim \text{Normal}(\mu_i, \sigma_i ),\;\; i=1,2$$

This is really why Bayesian analysis is great: we can specify what we think a fair price is through the $$\mu_i$$ parameter, and express uncertainty of our guess in the $$\sigma_i$$ parameter.

We'll assume two prizes per suite for brevity, but this can be extended to any number. The true price of the prize suite is then given by $$\text{Prize}_1 + \text{Prize}_2 + \epsilon$$, where $$\epsilon$$ is some error term.

We are interested in the updated $$\text{True Price}$$ given we have observed both prizes and have belief distributions about them. We can perform this using PyMC3.

Lets make some values concrete. Suppose there are two prizes in the observed prize suite:

1. A trip to wonderful Toronto, Canada!
2. A lovely new snowblower!

We have some guesses about the true prices of these objects, but we are also pretty uncertain about them. I can express this uncertainty through the parameters of the Normals:

\begin{align} & \text{snowblower} \sim \text{Normal}(3 000, 500 )\\\\ & \text{Toronto} \sim \text{Normal}(12 000, 3000 )\\\\ \end{align}

For example, I believe that the true price of the trip to Toronto is 12 000 dollars, and that there is a 68.2% chance the price falls 1 standard deviation away from this, i.e. my confidence is that there is a 68.2% chance the trip is in [9 000, 15 000].

In the first step we have to get samples from the posterior distribution.

My solution

I thought that we just have to sample values from the following distribution: $$p(x|{\rm obs}) \propto p({\rm obs}|x) p(x)$$ where the priror $$p(x)$$ is given by: $$p(x) = \frac{1}{\sqrt{2\pi}\tilde{\sigma} }\exp\left(-\frac{1}{2} \frac{(x-\tilde{\mu})^2}{\tilde{\sigma}^2} \right), \hspace{10.0mm} \tilde{\mu}=35000, \tilde{\sigma}=7500$$ The likelihood $$p({\rm obs}| x)$$ is obtained from the sum of normal distributions (snowblower $$\sim \mathcal{N}(\mu_1,\sigma_1)$$ and Toronto trip $$\sim \mathcal{N}(\mu_2,\sigma_2)$$):
$$p({\rm obs}|x) = \frac{1}{\sqrt{2\pi}\sigma_{12} } \exp\left(-\frac{1}{2} \frac{(x-\mu_{12})^2}{\sigma^2_{12}} \right), \hspace{10.0mm} \mu_{12}=\mu_1+\mu_2, \sigma_{12}=\sqrt{\sigma^2_1+\sigma^2_2 }\approx \sigma_2$$

By "completing the square" of the sum of arguments in the exponential functions in $$p(x|{\rm obs})$$, one can show that the maximum of the posterior is at: $${\rm argmax}_x \hspace{2mm} p(x|{\rm obs})= \frac{\mu_{12}\tilde{\sigma}^2 + \tilde{\mu}\sigma^2_{12} }{\sigma^2_{12}+\tilde{\sigma}^2 } \approx 17750$$

With the following code I obtain a $$p(x|{\rm obs})$$ with a peak at the predicted location

import pymc3 as pm
import numpy as np

mu_1, mu_2, mu_k = 3e3, 12e3, 35e3
mu_12= mu_1 + mu_2

sd_1, sd_2, sd_k = 0.5e3, 3e3, 7.5e3
sd_12= np.sqrt(sd_1**2 + sd_2**2) # approx equal to sd_2

with pm.Model() as model:

priror = pm.Normal("priror", mu=mu_k, sd=sd_k)
obs = pm.Normal("obs", mu=priror, sd=sd_12, observed=np.array([mu_12]))

### we can exchange the positions of priror and likelihood
### and get the same result
# priror = pm.Normal("priror", mu=mu_12, sd=sd_12)
# obs = pm.Normal("obs", mu=priror, sd=sd_k, observed=np.array([mu_k]))

trace = pm.sample(40000)


The problem is that in the textbook the posterior is peaked at $$x=20000$$. I would be glad if I get an idea from where the difference comes from.

Solution in the textbook

For comparison, the code in the textbook is given below:

import pymc3 as pm

data_mu = [3e3, 12e3]

data_std =  [5e2, 3e3]

mu_prior = 35e3
std_prior =  75e2
with pm.Model() as model:
true_price = pm.Normal("true_price", mu=mu_prior, sd=std_prior)

prize_1 = pm.Normal("first_prize", mu=data_mu[0], sd=data_std[0])
prize_2 = pm.Normal("second_prize", mu=data_mu[1], sd=data_std[1])
price_estimate = prize_1 + prize_2

logp = pm.Normal.dist(mu=price_estimate, sd=(3e3)).logp(true_price)
error = pm.Potential("error", logp)

trace = pm.sample(50000, step=pm.Metropolis())
burned_trace = trace[10000:]

price_trace = burned_trace["true_price"]