CODEWITHSUNDEEP
pymc3
Helmut Strey
Helmut Strey

Reputation: 100

How to create a container in pymc3

I am trying to build a model for the likelihood function of a particular outcome of a Langevin equation (Brownian particle in a harmonic potential):

Here is my model in pymc2 that seems to work: https://github.com/hstrey/BayesianAnalysis/blob/master/Langevin%20simulation.ipynb

#define the model/function to be fitted.
def model(x): 
    t = pm.Uniform('t', 0.1, 20, value=2.0)
    A = pm.Uniform('A', 0.1, 10, value=1.0)

    @pm.deterministic(plot=False)
    def S(t=t):
       return 1-np.exp(-4*delta_t/t)

    @pm.deterministic(plot=False)
    def s(t=t):
       return np.exp(-2*delta_t/t)

    path = np.empty(N, dtype=object)

    path[0]=pm.Normal('path_0',mu=0, tau=1/A, value=x[0], observed=True)
    for i in range(1,N):
        path[i] = pm.Normal('path_%i' % i,
                        mu=path[i-1]*s,
                        tau=1/A/S,
                        value=x[i],
                        observed=True)
        return locals()

mcmc = pm.MCMC( model(x) )
mcmc.sample( 20000, 2000, 10 )

The basic idea is that each point depends on the previous point in the chain (Markov chain). Btw, x is an array of data, N is its length, delta_t is the time step =0.01. Any idea how to implement this in pymc3? I tried:

# define the model/function for diffusion in a harmonic potential
DHP_model = pm.Model()
with DHP_model:
    t = pm.Uniform('t', 0.1, 20)
    A = pm.Uniform('A', 0.1, 10)

    S=1-pm.exp(-4*delta_t/t)

    s=pm.exp(-2*delta_t/t)

    path = np.empty(N, dtype=object)

    path[0]=pm.Normal('path_0',mu=0, tau=1/A, observed=x[0])
    for i in range(1,N):
        path[i] = pm.Normal('path_%i' % i,
                        mu=path[i-1]*s,
                        tau=1/A/S,
                        observed=x[i])

Unfortunately the model crashes as soon as I try to run it. I tried some pymc3 examples (tutorial) on my machine and this is working.

Thanks in advance. I am really hoping that the new samplers in pymc3 will help me with this model. I am trying to apply Bayesian methods to single-molecule experiments.

Upvotes: 2

Views: 750

Answers (2)

Noah
Noah

Reputation: 306

Rather than creating many individual normally-distributed 1-D variables in a loop, you can make a custom distribution (by extending Continuous) that knows the formula for computing the log likelihood of your entire path. You can bootstrap this likelihood formula off of the Normal likelihood formula that pymc3 already knows. See the built-in AR1 class for an example.

Since your particle follows the Markov property, your likelihood looks like

import theano.tensor as T

def logp(path):
    now = path[1:]
    prev = path[:-1]

    loglik_first = pm.Normal.dist(mu=0., tau=1./A).logp(path[0])
    loglik_rest = T.sum(pm.Normal.dist(mu=prev*ss, tau=1./A/S).logp(now))
    loglik_final = loglik_first + loglik_rest

    return loglik_final

I'm guessing that you want to draw a value for ss at every time step, in which case you should make sure to specify ss = pm.exp(..., shape=len(x)-1), so that prev*ss in the block above gets interpreted as element-wise multiplication.

Then you can just specify your observations with

path = MyLangevin('path', ..., observed=x)

This should run much faster.

Upvotes: 1

Helmut Strey
Helmut Strey

Reputation: 100

Since I did not see an answer to my question, let me answer it myself. I came up with the following solution:

# now lets model this data using pymc
# define the model/function for diffusion in a harmonic potential
DHP_model = pm.Model()
with DHP_model:
    D = pm.Gamma('D',mu=mu_D,sd=sd_D)
    A = pm.Gamma('A',mu=mu_A,sd=sd_A)

    S=1.0-pm.exp(-2.0*delta_t*D/A)

    ss=pm.exp(-delta_t*D/A)

    path=pm.Normal('path_0',mu=0.0, tau=1/A, observed=x[0])
    for i in range(1,N):
        path = pm.Normal('path_%i' % i,
                        mu=path*ss,
                        tau=1.0/A/S,
                        observed=x[i])

    start = pm.find_MAP()
    print(start)
    trace = pm.sample(100000, start=start)

unfortunately, this code takes at N=50 anywhere between 6hours to 2 days to compile. I am running on a pretty fast PC (24Gb RAM) running Ubuntu. I tried to using the GPU but that runs slightly slower. I suspect memory problems since it uses 99.8% of the memory when running. I tried the same calculation with Stan and it only takes 2min to run.

Upvotes: 1

Related Questions