# Help with pm.Potential

**URL:** <https://discourse.pymc.io/t/help-with-pm-potential/6562>\
**Category:** Questions\
**Created:** [January 11, 2021, 10:07am UTC](https://discourse.pymc.io/t/help-with-pm-potential/6562 "2021-01-11T10:07:37Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![sspickle](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/sspickle/32/3532_2.png) [@sspickle](https://discourse.pymc.io/u/sspickle)\
**Post date:** [January 11, 2021, 10:07am UTC](https://discourse.pymc.io/t/help-with-pm-potential/6562/1 "2021-01-11T10:07:37Z")

</div>

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](https://www.inference.org.uk/mackay/itila/book.html)" .

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:

```auto
    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

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [January 11, 2021, 6:16pm UTC](https://discourse.pymc.io/t/help-with-pm-potential/6562/2 "2021-01-11T18:16:26Z")

</div>

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`).

---

<div class="post-metadata">

**Author:** ![sspickle](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/sspickle/32/3532_2.png) [@sspickle](https://discourse.pymc.io/u/sspickle)\
**Post date:** [January 11, 2021, 10:37pm UTC](https://discourse.pymc.io/t/help-with-pm-potential/6562/3 "2021-01-11T22:37:36Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [January 12, 2021, 6:02am UTC](https://discourse.pymc.io/t/help-with-pm-potential/6562/4 "2021-01-12T06:02:44Z")

</div>

According to [Wikipedia](https://en.m.wikipedia.org/wiki/Truncated_distribution) 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:

```python

λ = 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))))

```

---

<div class="post-metadata">

**Author:** ![sspickle](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/sspickle/32/3532_2.png) [@sspickle](https://discourse.pymc.io/u/sspickle)\
**Post date:** [January 12, 2021, 9:34am UTC](https://discourse.pymc.io/t/help-with-pm-potential/6562/6 "2021-01-12T09:34:34Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [January 12, 2021, 9:56am UTC](https://discourse.pymc.io/t/help-with-pm-potential/6562/7 "2021-01-12T09:56:00Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [January 12, 2021, 10:04am UTC](https://discourse.pymc.io/t/help-with-pm-potential/6562/8 "2021-01-12T10:04:26Z")

</div>

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](https://mc-stan.org/docs/2_20/reference-manual/sampling-statements-section.html#truncated-distributions)

---

<div class="post-metadata">

**Author:** ![sspickle](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/sspickle/32/3532_2.png) [@sspickle](https://discourse.pymc.io/u/sspickle)\
**Post date:** [January 12, 2021, 10:40am UTC](https://discourse.pymc.io/t/help-with-pm-potential/6562/9 "2021-01-12T10:40:55Z")

</div>

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

```auto
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:

```auto
    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?

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [January 12, 2021, 10:58am UTC](https://discourse.pymc.io/t/help-with-pm-potential/6562/10 "2021-01-12T10:58:30Z")

</div>

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

Edit: There is one example with censored data in this Notebook ([resources/DataAnalysis.ipynb at master · pymc-devs/resources · GitHub](https://github.com/pymc-devs/resources/blob/master/BCM/ParameterEstimation/DataAnalysis.ipynb)), 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](https://colab.research.google.com/drive/1tPcVA_xBtGEIH_eIxCNJC51yKQKZt3Tx?usp=sharing) (?)

---

<div class="post-metadata">

**Author:** ![sspickle](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/sspickle/32/3532_2.png) [@sspickle](https://discourse.pymc.io/u/sspickle)\
**Post date:** [January 13, 2021, 2:13am UTC](https://discourse.pymc.io/t/help-with-pm-potential/6562/11 "2021-01-13T02:13:48Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![sspickle](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/sspickle/32/3532_2.png) [@sspickle](https://discourse.pymc.io/u/sspickle)\
**Post date:** [January 13, 2021, 2:16am UTC](https://discourse.pymc.io/t/help-with-pm-potential/6562/12 "2021-01-13T02:16:55Z")

</div>

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

```auto
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)

```

```auto
	mean	sd hdi_3%	hdi_97%	
λ	2.895	0.159	2.603	3.199

```

Thank you!