CODEWITHSUNDEEP
pythonpymcmcmc
r_31415
r_31415

Reputation: 8982

How should I use @pm.stochastic in PyMC?

Fairly simple question: How should I use @pm.stochastic? I have read some blog posts that claim @pm.stochasticexpects a negative log value:

@pm.stochastic(observed=True)
def loglike(value=data):
  # some calculations that generate a numeric result
  return -np.log(result)

I tried this recently but found really bad results. Since I also noticed that some people used np.log instead of -np.log, I give it a try and worked much better. What is really expecting @pm.stochastic? I'm guessing there was a small confusion on the sign required due to a very popular example using something like np.log(1/(1+t_1-t_0)) which was written as -np.log(1+t_1-t_0)

Another question: What is this decorator doing with the value argument? As I understand it, we start with some proposed value for the priors that need to enter in the likelihood and the idea of @pm.stochastic is basically produce some number to compare this likelihood to the number generated by the previous iteration in the sampling process. The likelihood should receive the value argument and some values for the priors, but I'm not sure if this is all value is doing because that's the only required argument and yet I can write:

@pm.stochastic(observed=True)
def loglike(value=[1]):
  data = [3,5,1] # some data
  # some calculations that generate a numeric result
  return np.log(result)

And as far as I can tell, that produces the same result as before. Maybe, it works in this way because I added observed=True to the decorator. If I would have tried this in a stochastic variable with observed=False by default, value would be changed in each iteration trying to obtain a better likelihood.

Upvotes: 1

Views: 1175

Answers (2)

RichardBJ
RichardBJ

Reputation: 377

Yes confusion is easy since the @stochastic returns a likelihood which is the opposite of the error essentially. So you take the negative log of your custom error function and return THAT as your log-likelihood.

Upvotes: 0

Abraham D Flaxman
Abraham D Flaxman

Reputation: 2979

@pm.stochastic is a decorator, so it is expecting a function. The simplest way to use it is to give it a function that includes value as one of its arguments, and returns a log-likelihood.

You should use the @pm.stochastic decorator to define a custom prior for a parameter in your model. You should use the @pm.observed decorator to define a custom likelihood for data. Both of these decorators will create a pm.Stochastic object, which takes its name from the function it decorates, and has all the familiar methods and attributes (here is a nice article on Python decorators).

Examples:

A parameter a that has a triangular distribution a priori:

@pm.stochastic
def a(value=.5):
    if 0 <= value < 1:
        return np.log(1.-value)
    else:
        return -np.inf

Here value=.5 is used as the initial value of the parameter, and changing it to value=1 raises an exception, because it is outside of the support of the distribution.

A likelihood b that has is normal distribution centered at a, with a fixed precision:

@pm.observed
def b(value=[.2,.3], mu=a):
    return pm.normal_like(value, mu, 100.)

Here value=[.2,.3] is used to represent the observed data.

I've put this together in a notebook that shows it all in action here.

Upvotes: 3

Related Questions