Why am I getting Bad initial energy with this simple model

Questions
1

Hi,

I’m testing a simple model, where the prior for theta is Uniform in 0,1 and the likelihood is uniform in 0 theta. When I run the code, I obtain SamplingError("Bad initial energy").

with pm.Model() as basic_model:
    theta = pm.Uniform('theta', lower=0, upper=1)

    y = pm.Uniform('y', lower=0, upper=theta,
                   observed=[0.49131995252993826, 0.2774397121516236, 0.5381935236905224, 0.19753121107715765])

    print(basic_model.check_test_point())

    trace = pm.sample(1000, progressbar=True, chains=1, random_seed=1)

The output for check_test_point is

theta_interval__ -1.39
y -inf
Name: Log-probability of test_point, dtype: float64

What am I doing wrong?

Thanks in advance

2

Hi,

I can run your model, despite the fact that it seems pretty odd, namely to use a Uniform distribution as the Likelihood. Due to that fact, I get a lot of divergences.

3

Hi @luisroque

Please, can you run it like this?

with pm.Model() as basic_model:
    theta = pm.Uniform('theta', lower=0, upper=1)

    y = pm.Uniform('y', lower=0, upper=theta,
                   observed=[0.49131995252993826, 0.2774397121516236, 0.5381935236905224, 0.19753121107715765])

    print(basic_model.check_test_point())

    trace = pm.sample(1000, progressbar=True, chains=1, random_seed=1)

With that dataset, check_test_point returns the following:

theta_interval__ -1.39
y -inf
Name: Log-probability of test_point, dtype: float64

and when I do the sampling I get SamplingError("Bad initial energy").

I´ll update the original questions to use this code.

By the way, I took the model from the following example in the book Introduction to Probability by Bertsekas and Tsitsiklis (pp 414)

Romeo and Juliet start dating, but Juliet will be late on any
date by a random amount X, uniformly distributed over the interval [0,
theta]. The parameter theta is unknown and is modeled as the value of a random variable,
uniformly distributed between zero and one hour. Assuming that Juliet was late by
an amount x on their first date, how should Romeo use this information to update
the distribution of theta?

I can solve this problem analytically and I want to compare the analytical result with the one obtained by PyMC3.

Thanks

4

Hm, I see. It still looks strange to me. But nevertheless, from what I can see, theta has to be bigger than the biggest value from your data. If it isn’t the log-prob goes to infinity. You can also see that by plotting the posterior for theta. The problem is that PyMC3 does not know that and the initial value that it sets up ends up being smaller. You can define the initial point yourself using testval.

with pm.Model() as basic_model:
    theta = pm.Uniform('theta', lower=0, upper=1, testval=0.54)

    y = pm.Uniform('y', lower=0, upper=theta,
                   observed=[0.49131995252993826, 0.2774397121516236, 0.5381935236905224, 0.19753121107715765])

    print(basic_model.check_test_point())

    trace = pm.sample(1000, progressbar=True, chains=1, random_seed=1)

See the domain for the posterior of theta:

image

5

I think @luisroque is right about the testval being needed to ensure a valid starting energy.

In this kind of models it is (unfortunately) usually necessary to set a prior that is already compatible with the observed data. In this case you know that theta cannot be lower than your highest observation, so you can try something like

with pm.Model() as basic_model:
    theta = pm.Uniform('theta', lower=0.5381935236905224, upper=1)

    y = pm.Uniform('y', lower=0, upper=theta,
                   observed=[0.49131995252993826, 0.2774397121516236, 0.5381935236905224, 0.19753121107715765])
6

Ok, thanks. I’ve modified the code to do that and now it works fine.

data = [0.49131995252993826, 0.2774397121516236, 0.5381935236905224, 0.19753121107715765]

with pm.Model() as basic_model:
    theta = pm.Uniform('theta', lower=max(data), upper=1)

    y = pm.Uniform('y', lower=0, upper=theta,
                   observed=data)

    print(basic_model.check_test_point())

    trace = pm.sample(1000, progressbar=True)
1 Like