0

If you have a random variable $X$ and a function $f$, you can define $y=f(X)$ as a new random variable with a probability density function as follows:

$p(y)=(f^{-1})'(y)p(x)$. For details see here.

Now I have defined a random variable alpha, with an exponential distribution in the following code. I want to add to my model, log(alpha) as a new random variable. How should I implement it in my model?

I already made an effort but it seems that it is wrong, and the reason as been pointed out in answers is the fact that I used stochastic decorator rather than a deterministic one. But because later I want to apply MCMC Metropolis on this variable I need it to be statistic! To clarify it more, I want to apply a Gaussian proposal on the log(alpha). So I need to hand in an stochastic input to Metropolis function.

So this is my model:

import numpy as np
import pymc
lambd=1;
__all__=['alpha']
alpha=pymc.Exponential('alpha', beta=lambd)

@pymc.stochastic(plot=False)
def logalpha(value=0,c=alpha):
    return np.log(c)
Cupitor
  • 11,007
  • 19
  • 65
  • 91
  • Is `logalpha` supposed to the the `log` of the stochastic variable `alpha` or are you intending it to be an independent variable? – bogatron Oct 17 '13 at 13:49
  • I want it to represent the log of the stochastic variable alpha, but I want to be able to sample it as well. – Cupitor Oct 17 '13 at 14:03

2 Answers2

1

log alpha is a deterministic function of your alpha, so you should model it as @deterministic. A good toy example that mirrors your own problem is the regression example.

Max
  • 1,670
  • 1
  • 12
  • 17
  • the problem is when I make it deterministic. I cannot sample from it. – Cupitor Oct 17 '13 at 14:01
  • PyMC should accrue samples even for deterministic variables. In your example, `model.trace('logalpha')[:]` should give you the samples. – Max Oct 17 '13 at 14:37
  • I edited my question. What I need is to be able to apply MCMC metropolis on log(alpha) not on alpha! – Cupitor Oct 17 '13 at 16:03
  • @Naji, sampling from alpha, and applying log() is the same as sampling directly from log(alpha) – Cam.Davidson.Pilon Oct 17 '13 at 16:32
  • @Cam.Davidson.Pilon, I do understand that. Because I am using a Gaussian proposal and I am applying it on log(alpha), I need to apply my sampling method on log(alpha) using "use_step_method" – Cupitor Oct 17 '13 at 16:42
  • What is your ultimate goal? Why do you need to use Gaussian proposals? – Cam.Davidson.Pilon Oct 17 '13 at 16:46
  • @Cam.Davidson.Pilon, Because exponential function is unbounded, however I am trying to use a Gaussian proposal to estimate it by metropolis. I don't want the PyMC to sample it for me. – Cupitor Oct 17 '13 at 17:28
0

As @Max already mentioned, logalpha should be a deterministic variable, since it's value is uniquely determined by alpha. Whenever your model is sampled, the value of logalpha will be updated accordingly. For example:

>>> import numpy as np
>>> import pymc
>>> lambd = 1
>>> 
>>> alpha = pymc.Exponential('alpha', beta=lambd)
>>> 
>>> @pymc.deterministic(plot=False)
... def logalpha(value=0, c=alpha):
...     return np.log(c)
... 
>>> M = pymc.Model([alpha, logalpha])
>>> for i in range(3):
...     M.draw_from_prior()
...     print (alpha.value, logalpha.value)
... 
(array(1.888410537018971), 0.63573548954043602)
(array(0.23180935966225977), -1.4618399707110767)
(array(0.3381518219555991), -1.0842603069656513)
bogatron
  • 18,639
  • 6
  • 53
  • 47
  • thank you for your reply, but I want to be able to sample log(alpha) itself using MCMC Metropolis with a custom proposal distribution(in my case Gaussian) – Cupitor Oct 17 '13 at 16:03
  • I don't see how this prevents you from doing that. But you could also just create your variable as a pymc `lognormal_like` variable (http://pymc-devs.github.io/pymc/distributions.html#continuous-distributions). – bogatron Oct 17 '13 at 16:20
  • A deterministic variable cannot be handed as a parameter to use_step_method method. I think you cannot define a random "variable" from a lognormal_like "distribution". Even if you could that would be different from the alpha random variable I defined. – Cupitor Oct 17 '13 at 16:22