CODEWITHSUNDEEP
pythonpymc
Chris L. Barnes
Chris L. Barnes

Reputation: 482

Putting bounds on stochastic variables in PyMC

I have a variable A which is Bernoulli distributed, A = pymc.Bernoulli('A', p_A), but I don't have a hard value for p_A and want to sample for it. I do know that it should be small, so I want to use an exponential distribution p_A = pymc.Exponential('p_A', 10).

However, the exponential distribution can return values higher than 1, which would throw off A. Is there a way of bounding the output of p_A without having to re-implement either the Bernoulli or the Exponential distributions in my own @pymc.stochastic-decorated function?

Upvotes: 1

Views: 868

Answers (3)

Maria
Maria

Reputation: 108

PyMC provides bounds. The following should also work:

p_A = pymc.Bound(pymc.Exponential, upper=1)('p_A', lam=10)

Upvotes: 1

Chris L. Barnes
Chris L. Barnes

Reputation: 482

For any other lost souls who come across this:

I think the best solution for my purposes (that is, I was only using the exponential distribution because the probabilities I was looking to generate were probably small, rather than out of mathematical convenience) was to use a Beta function instead.

For certain parameter values it approximates the shape of an exponential function (and can do the same for binomials and normals), but is bounded to [0 1]. Probably only useful for doing things numerically, though, as I imagine it's a pain to do any analysis with.

Upvotes: 0

Ramon Crehuet
Ramon Crehuet

Reputation: 4017

You can use a deterministic function to truncate the Exponential distribution. Personally I believe it would be better if you use a distribution that is bound between 0 and 1, but to exactly solve your problem you can do as follows:

import pymc as pm
p_A = pm.Exponential('p_A',10)

@pm.deterministic
def p_B(p=p_A):
    return min(1, p)

A = pm.Bernoulli('A', p_B)

model = dict(p_A=p_A, p_B=p_B, A=A)
S = pm.MCMC(model)
S.sample(1000)
p_B_trace = S.trace('p_B')[:]

Upvotes: 3

Related Questions