Why can't PyMC3 fit a uniform distribution with a Normal prior?

Question

My question might be slightly ill-posed. My confusion comes from the oft repeated "Markov Chain Monte Carlo is a technique to solve the problem of sampling from a complicated distribution." -- http://jeremykun.com/2015/04/06/markov-chain-monte-carlo-without-all-the-bullshit/

So I tried to see if I could learn a stupidly simple distribution (uniform). I created some fake data drawn uniformly from (0-5)

observed_data = np.random.uniform(low=0, high=5, size=1000)

that looks reasonable: .

My plan was to use PyMC3 to fit this distribution -- but starting with a Normal distribution. I know you're thinking hold up, that isn't right, but I was under the impression that a Normal distribution would just be the prior that MCMC would be flexible enough to discover the underlying distribution. So, doing that:

with pm.Model():
    mu = pm.Normal('mu', 0, 1)
    sigma = 1.
    returns = pm.Normal('returns', mu=mu, sd=sigma, observed=observed_data)

    step = pm.NUTS()
    trace = pm.sample(15000, step)

sns.distplot(trace[-5000:]['mu'], label='PyMC3');
plt.legend()
plt.savefig("sampled.png")

And the resulting trace: . Looking at this plot I can see that PyMC3 has fit a Normal distribution (has it?). It looks to me that it has decided that it's pretty confident that the center of this distribution is ~2.45 +/- .1. This surprises me!

I want to know two things:

1) Where is my conceptual blind-spot here? Why am I surprised at this result?

2) Is there a way for getting a distribution that looks similar to the observed data (the uniform distribution in the first plot)? How?

Edit/addition

To model the actual predictive distribution, I made two changes to the above code. Made sigma a HalfNormal distribution and made a predictive call.

with pm.Model():
    mu = pm.Normal('mu', 0, 1)
    # sigma = 1.0
    sigma = pm.HalfNormal('sigma', sd=2.0)

    predictive = pm.Normal('predictive', mu=mu, sd=sigma)
    returns = pm.Normal('returns', mu=mu, sd=sigma, observed=observed_data)

    step = pm.NUTS()
    trace = pm.sample(500000, step, tune=20000)

Looking at the observed data (uniform distribution) compared with the predictive distribution (Normal), gives the following:

.

As Björn suggests, the Normal isn't flexible enough to mimic the uniform (though it does look like it's being deformed in that way).

Answer 1

I think the program has done exactly what you asked it to and has done so pretty well (from the limited information you show). Given the information that you are 100% certain that this is data from a normal distribution with standard deviation 1 and with a mean that you are pretty sure is not that far away from 1, it has identified the posterior distribution for the mean of this normal distribution. Unsurprisingly this is close to the mean of the distribution you simulated from.

The overall distribution it thinks it saw would actually be obtained by simulating from a normal with this distribution of means and a standard deviation of 1 (that would be the posterior predictive distriution describing how a new observation from the same distribution is predicted to be distributed like). That would of course still be a normal distribution (i.e. with a peak around the mean) and dropping off towards the edge of the uniform distribution (because you told it that sigma=1). No matter what you do a single normal distribution is too unflexible to approximate a uniform distribution all that well. Perhaps a mixture of normals (with unknown means and standard deviations) would look a bit better, but if you told the program that you have drawn the data from a uniform distribution, then I would expect the best result.

Answer 2

I can't really read the code, but using a normal model to approximate a uniform should always result in exactly this result: a normal centered at the uniform mean with variance to match. If you want a better approximation, give the model more flexibility to adapt to the data. If you want the best outcome, though, you'll have to correctly specify the model.

I think that the root of your understanding is in parsing the quotation

My confusion comes from the oft repeated "Markov Chain Monte Carlo is a technique to solve the problem of sampling from a complicated distribution."

Note the the quote says "complicated," not "unknown." In Bayesian statistics, we're perfectly capable of writing out what the probability distribution over the parameters given the data must be -- but it's usually completely intractable to evaluate. Typically, the math is only straightforward for the simplest models, and advanced calculus can take us a little further, but a Bayesian problem of even modest complexity requires simulation to resolve, particularly Bayesian hierarchical models.

Your writing here also supports my interpretation of where the misunderstanding lies.

So I tried to see if I could learn a stupidly simple distribution (uniform). I created some fake data drawn uniformly from (0-5).

Again, I think that you've miscast the problem. PyMC isn't "learning" a distribution. You're just drawing samples from a specified functional form. The way of going about that is somewhat complicated, but it's just about drawing samples from the model that you've specified. The normal has support over the real line, so PyMC is drawing samples from the reals. The data more strongly support samples drawn from $[0,5]$ (can you see why?) but nothing about the model you've specified excludes $\pm6$ (for example) from having positive probability.

I would recommend reading more about Bayesian statistics. Gelman's Bayesian Data Analysis is a great place to start.