test.py

import pymc3 as pm
import numpy as np
from pymc3.math import switch
import theano
from matplotlib import pyplot as plt
import pandas as pd
from scipy import optimize
from interpolate_marketdata import *
import theano.tensor as T
N_sims = 1000

data = [13,24,8,24,7,35,14,11,15,11,22,22,11,57,
11,19,29,6,19,12,22,12,18,72,32,9,7,13,19,
23,27,20,6,17,13,10,14,6,16,15,7,2,15,15,19,
70,49,7,53,22,21,31,19,11,18,20,12,35,17,23,
17,4,2,31,30,13,27,0,39,37,5,14,13,22]

""" detect any regime change in factors 2 & 3 """
n_count_data = len(data)
basic_model = pm.Model()

with basic_model:
  mu_1 = pm.Normal('mu_1', mu=0, sd=10, shape=1)
  mu_2 = pm.Normal('mu_2', mu=0, sd=10, shape=1)
  sigma_1 = pm.HalfNormal('sigma_1', sd=1, shape=1)
  sigma_2 = pm.HalfNormal('sigma_2', sd=1, shape=1)
  tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data)

  mu1_print = T.printing.Print('mu_1')(mu_1)
  mu2_print = T.printing.Print('mu_2')(mu_2)
#   sd_print = T.printing.Print('sd')(sd)

  idx = np.arange(n_count_data)
  mu_ = switch(tau >= idx, mu_1, mu_2)
  sigma_ = switch(tau >= idx, sigma_1, sigma_2)
  # observed data can be a ndarray or a pandas dataframe
  obs = pm.Normal('obs', mu=mu_, sd=sigma_, observed=data)
  # obtain starting values via MAP
  start = pm.find_MAP(fmin=optimize.fmin_powell)
  step=pm.NUTS()
  trace = pm.sample(N_sims,step,progressbar=True,start=start)
""" get samples from posterior distribution """
pm.traceplot(trace)
plt.show()