How to use free random variables in step method function?

Questions
1

Hi,

I am trying to write a code where I pass a free random variable to a custom step method in order to calculate the new proposal value (in other words, the new proposal value is depending on the current value of the free random variable). I am doing a single site update, meaning that in my array q of length 100, I chose a position randomly and then only update only that position at a time.

I have tried various ways, e.g.
ni = some constant
R = Array of length 100

with Model() as model:
    gamma = pm.Gamma('gamma', 81, 100)
    I = MyDistr('I', gamma=gamma, ni=ni, testval=0.5, shape=100)
    updI = MyStepMethod(vars=[I], gamma=gamma, ni=ni)
    trace = pm.sample(10000, step=[updI], tune=2000)

class MyStepMethod(ArrayStep):
def __init__(self, vars, R, ni, gamma, model=None):
    self.vars = vars
    self.R = R
    self.ni = ni
    self.gamma = gamma
    model = pm.modelcontext(model)
    super(MyStepMethod, self).__init__(vars, [model.fastlogp])

def astep(self, q0, logp):

    q = q0.copy()
    ind = np.random.randint(low=0, high=self.ni, size=1)
    q[ind] = self.R[ind] - self.gamma  # setting new value
    lognew = logp(q)
    logalt = logp(q0)   
    # an then do a Metroplois Hastings with logold and lognew

The error message I get is: ValueError: setting an array element with a sequence. I dont understand this, because I thought gamma is not a sequence but rather a 1 dim element…?

I also tried passing gamma as vars parameter to the stepmethod function:

updI = MyStepMethod(vars=[I, gamma], R=R, ni=ni)

and then the idea is to use the last element of q as my gamma value:

def astep(self, q0, logp):

    q = q0.copy()
    ind = np.random.randint(low=0, high=self.ni, size=1)
    q[ind] = self.R[ind] - q[-1]  # setting new value
    lognew = logp(q[:-1])
    logalt = logp(q0[:-1])

However, I get the error message already at the updI = MyStepMethod(vars=[I, gamma] … line saying ValueError: Shape of variable not known gamma.
Also, setting gamma = pm.Gamma(‘gamma’, 81, 100, shape=1) does not resolve the problem…

Does anyone has any suggestion on how to go about it? Thank you!

2

Did you try with a distribution from pymc3? The error could also related to your implementation of MyDistr.

3

Yes, I have tried it with a pymc3 distribution. I came up with an easy example (similar to this ) where it actually works with a Normal distribution as prior for the mean parameter.

import pymc3 as pm
import numpy as np
import pymc3.backends.ndarray
import matplotlib.pyplot as plt
from pymc3.step_methods.arraystep import ArrayStep
from pymc3 import Continuous

class MyStepMethod(ArrayStep):

    def __init__(self, vars, model=None):
        self.vars = vars
        model = pm.modelcontext(model)
        super(MyStepMethod, self).__init__(vars, [model.fastlogp])
        self.accrate = 0
        self.count = 0

    def astep(self, q0, logp):

        q = q0.copy()
        q[0] = np.random.normal(q[-1], 1) # using current gamma value as mean

        logalt, logneu = logp(q0), logp(q)
        accept = logneu - logalt

        if (accept > 0) or (np.random.uniform(low=0, high=1) < np.exp(accept)):
            self.accrate += 1
            return q

        else: return q0

with pm.Model() as model:
    mean = pm.Normal("test", 5, 0.001)
    I = pm.Normal("I", mean,1)

    updI = MyStepMethod(vars=[I, mean])
    step1 = pm.NUTS(vars=[mean])
    trace = pm.sample(5000, [updI, step1], start={'I':0.5}, njobs=1, chains=1, tune=1000)

pm.traceplot(trace)
plt.show()

However, if I use a different distribution for my gamma parameter and e.g. change the line

mean = pm.Normal(“test”, 5, 0.001)

to

mean = pm.Gamma(“test”, 81, 100)

I get the shape error:

Traceback (most recent call last):
  File "/home/test.py", line 50, in <module>
    updI = MyStepMethod(vars=[I, mean])
  File "/home/test.py", line 26, in __init__
    super(MyStepMethod, self).__init__(vars, [model.fastlogp])
  File "/home/user/.local/lib/python3.6/site-packages/pymc3/step_methods/arraystep.py", line 113, in __init__
    self.ordering = ArrayOrdering(vars)
  File "/home/user/.local/lib/python3.6/site-packages/pymc3/blocking.py", line 36, in __init__
    raise ValueError('Shape of variable not known %s' % name)
ValueError: Shape of variable not known test

Any idea on why it only works for normal distributions or how I can get it to work for a gamma distributed variable?

4

In your code you are doing:

    updI = MyStepMethod(vars=[I, mean])
    step1 = pm.NUTS(vars=[mean])

As you are updating RV mean twice. I am not sure that this is a valid thing to do - is ergodicity still satisfied?

Otherwise, the reason that you are getting an error for Gamma but not for normal or other distribution defined on (-\infty, \infty) is because PyMC3 automatically transformed all bounded parameters to unbounded internally, thus the free parameter the sampler are actually “seeing” is the transformed parameter, in this case would be test_log__, which is the log transformed of test. A quick fix in this case is to do

with pm.Model() as model:
    mean = pm.Gamma("test", 81, 100, transform=None)
    I = pm.Normal("I", mean,1)

which turn off the auto transformation.

5

Thanks very much for you answer!

I agree, the RV mean is updated twice, but since I am not changing the value of mean in MyStepMethod, ergodicity should still be satisfied. Though, I wonder if there is a different way of passing the current value of mean into the MyStepMethod without using it as an update step, but rather to use it as a variable for calculating the new proposal value for I?

6

Oh I see. Yes there are ways to pass the current value, the trick is to compile a new logp function that conditioned on the rest of the value that the current sampler is not interested in, using the model.logp_dlogp_function:

You can have a look at


as example