Help with pm.Potential

Questions
1

I’m trying to use PyMC3 to work out an example (exercise 3.3, p 47) from David MacKay’s great book " Information Theory, Inference, and Learning Algorithms" .

I was able to use a crude (manual) grid search approach and it worked OK. I was hoping to learn more about PyMC3 by using it for the same purpose. Here’s what I tried:

import numpy as np
import pymc3 as pm
import arviz as az
import theano.tensor as tt

print(f"Running on PyMC3 v{pm.__version__}")

np.random.seed(451)
x = np.random.exponential(3,size=500)

minx=1
maxx=20

obs = x[np.where(~((x<minx) | (x>maxx)))] # remove values outside range

with pm.Model() as expModel:
    λ = pm.Exponential("λ", lam=5)
    x = pm.Exponential('x', lam=1/λ, observed=obs)
    clip_exp = pm.Potential("clip_exp", tt.switch(((x<minx)|(x>maxx)),-np.inf,0))
    trace= pm.sample(2000, tune=1000, return_inferencedata=True)

az.summary(trace)

Sadly it doesn’t appear to have worked:

    mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_mean  ess_sd 
λ  3.683  0.188   3.333    4.051      0.003    0.002    3268.0  3268.0

Any suggestions would be appreciated!

thank you,
-Steve

1 Like
2

I think you want to use a bound on lambda (or 1/lambda)?

x is not a random variable since you observed it, so the pm.Potential is not doing anything (x will always equal obs, and the clip_exp will always return 0).

3

Ah, OK. Thank you for that insight. The idea is that I’m observing x, which is exponentially distributed, but I’m only able to observe values of x between minx, and maxx. Values outside this range are not available to me. In my model \lambda has the same units as x, so I need to specific lam (which has units of 1/x) as the reciprocal of \lambda. I looked into pm.Bound, but my understanding is that it cannot be used for observed values. Is there a simple way to describe observed values that are censored in this way?

thanks,
-Steve

1 Like
4

According to Wikipedia it seems you need to add a correction term equal to cdf(max) - cdf(min). My guess would be a potential like this might do the job:


λ = pm.Exponential("λ", lam=5)
x = pm.Exponential('x', lam=1/λ, observed=obs)    

exp_dist = pm.Exponential.dist(lam=1/λ). # this is not part of the model, just used to get the logcdf
clip_exp = pm.Potential("clip_exp", tt.log(tt.exp(exp_dist.logcdf(maxx)) - tt.exp(exp_dist.logcdf(minx))))
1 Like
6

Ah, I see. It’s not independent of \lambda, just x. This didn’t seem to work, but I’ll poke around a bit more. Thanks for the help.

7

I wrongly used the argument mu in exp_dist instead of lam in my code snippet above (changed it now). Maybe that’s why it was not working?

8

The other thing is whether the value in the Potential should be negated (I never remember this).

This STAN page may be helpful to understand the general approach: 7.4 Sampling Statements | Stan Reference Manual

9

Thanks. Here’s what I’ve got right now:

import numpy as np
import pymc3 as pm
import arviz as az
import theano.tensor as tt

print(f"Running on PyMC3 v{pm.__version__}")

np.random.seed(451)
x = np.random.exponential(3,size=500)

minx=1
maxx=20

obs = x[np.where(~((x<minx) | (x>maxx)))] # remove values outside range

with pm.Model() as expModel:
    λ = pm.Exponential("λ", lam=1/5)  # prior exponential with mean of 5
    x = pm.Exponential('x', lam=1/λ, observed=obs) # obs exponential with mean of $\lambda$.
    
    exp_dist = pm.Exponential.dist(lam=1/λ) # this is not part of the model, just used to get the logcdf
    clip_exp = pm.Potential("clip_exp", -pm.math.logdiffexp(exp_dist.logcdf(maxx), exp_dist.logcdf(minx)))

    trace= pm.sample(2000, tune=1000, return_inferencedata=True)

az.summary(trace)

Still getting:

    mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_mean  ess_sd 
λ  3.874  0.206   3.484    4.257      0.004    0.003    3279.0  3279.0

I’m expecting \lambda to be around 3.

The -pm.math... doesn’t seem to matter (+ or -, makes no difference). Not sure how to approach debugging this. Any ideas?

10

Maybe the Potential needs to be multiplied by the number of observations (I am running out of ideas here :stuck_out_tongue: )?

Edit: There is one example with censored data in this Notebook (resources/DataAnalysis.ipynb at master · pymc-devs/resources · GitHub), section 5.5. It strengthens my conviction that the the correction term has to be proportional to the number of observations.

clip_exp = pm.Potential("clip_exp", pm.math.logdiffexp(exp_dist.logcdf(maxx), exp_dist.logcdf(minx)) * x.size)

Edit: It seems to work be working in this Colab NoteBook (?)

1 Like
11

Oh, wow. I never would’ve guessed that! Thank you. I’ll give that a try.

12

I think that worked! Here’s what I have now:

import numpy as np
import pymc3 as pm
import arviz as az
import theano.tensor as tt

print(f"Running on PyMC3 v{pm.__version__}")

np.random.seed(451)
x = np.random.exponential(3,size=500)

minx=1
maxx=20

obs = x[np.where(~((x<minx) | (x>maxx)))] # remove values outside range

def censor_correct(λ):
    result = -tt.log(tt.exp(-minx/λ)- tt.exp(-maxx/λ))*len(obs)
    return result
    
with pm.Model() as expModel:
    λ = pm.Exponential("λ", lam=1/5)  # prior exponential with mean of 5
    clip_exp = pm.Potential("clip_exp", censor_correct(λ))
    x = pm.Exponential('x', lam=1/λ, observed=obs) # obs exponential with mean of $\lambda$.
    
    trace= pm.sample(2000, tune=1000, return_inferencedata=True)

az.summary(trace)

	mean	sd	    hdi_3%	hdi_97%	
λ	2.895	0.159	2.603	3.199

Thank you!

4 Likes