CODEWITHSUNDEEP
pymc
spencerlyon2
spencerlyon2

Reputation: 9676

How to implement custom Gibbs sampling scheme in pymc

I have a hidden Markov stochastic volatility model (represented as a linear state space model). I am using a hand-written Gibbs sampling scheme to estimate parameters for the model. The actual sampler requires some fairly sophisticated update rules that I believe I need to write by hand. You can see an example of a Julia version of these update rules here.

My question is the following: how can I specify the model in a custom way and then hand the job of running the sampler and collecting the samples to pymc? In other words, I am happy to provide code to do all the heavy lifting (how to update each block of parameters on each scan -- utilizing full conditionals within each block), but I want to let pymc handle the "accounting" for me.

I realize that I will probably need to provide more information so that others can answer this question. The problem is I am not sure exactly what information will be useful. So, if you feel you can help me out with this, but need more information -- please let me know in a comment and I will update the question.

Upvotes: 2

Views: 1006

Answers (1)

Abraham D Flaxman
Abraham D Flaxman

Reputation: 2979

Here is an example of a custom sampler in PyMC2:

class BDSTMetropolis(mc.Metropolis):
    def __init__(self, stochastic):
        mc.Metropolis.__init__(self, stochastic, scale=1., proposal_sd='custom',
            proposal_distribution='custom', verbose=None, tally=False)

    def propose(self):
        T = self.stochastic.value

        T.u_new, T.v_new = T.edges()[0]
        while T.has_edge(T.u_new, T.v_new):
            T.u_new, T.v_new = random.choice(T.base_graph.edges())

        T.path = nx.shortest_path(T, T.u_new, T.v_new)
        i = random.randrange(len(T.path)-1)
        T.u_old, T.v_old = T.path[i], T.path[i+1]

        T.remove_edge(T.u_old, T.v_old)
        T.add_edge(T.u_new, T.v_new)
        self.stochastic.value = T

    def reject(self):
        T = self.stochastic.value
        T.add_edge(T.u_old, T.v_old)
        T.remove_edge(T.u_new, T.v_new)
        self.stochastic.value = T

It pretty different than your model, but it should demonstrate all the parts. Does that give you enough to go on?

Upvotes: 1

Related Questions