How to use the pymc.MCMC function in pymc

def fit_model(model):
    """ Fit model with MCMC, starting from MAP as initial value, and
    plot results"""

    #print "Fitting vars:"
    #print model.vars.describe()

    import time
    start_time = time.time()
    model.map = mc.MAP(model.vars)
    model.map.fit(method='fmin_powell', verbose=0)

    model.mcmc = mc.MCMC(model.vars)
    model.mcmc.use_step_method(mc.AdaptiveMetropolis, model.vars['gamma'])
    model.mcmc.sample(20000, 10000, 100, verbose=False, progress_bar=False)
    model.mcmc.wall_time = time.time() - start_time

sample_number: number of sampling. It must be greater than burn, however there is no check.
    burn: number of samples being burned.
    count1: observed counts of condition 1
    count2: observed counts of condition 2

    return: list of log2-ratios
    """
    lam1 = pymc.Uniform('U1',0,10000)   # prior of lambda is uniform distribution
    lam2 = pymc.Uniform('U2',0,10000)   # prior of lambda is uniform distribution    
    poi1 = pymc.Poisson('P1',lam1,value=count1,observed=True) # Poisson with observed value count1
    poi2 = pymc.Poisson('P2',lam2,value=count2,observed=True) # Poisson with observed value count2
    @deterministic
    def ratio (l1=lam1,l2=lam2):
        return log(l1,2) - log(l2,2)
    mcmcmodel  = pymc.MCMC([ratio,lam1,poi1,lam2,poi2])
    mcmcmodel.use_step_method(pymc.AdaptiveMetropolis,[ratio,lam1,lam2,poi1,poi2], delay=20000)
    if PROGRESS_BAR_ENABLED:
        mcmcmodel.sample(iter=sample_number, progress_bar=False, burn=burn)    
    else:
        mcmcmodel.sample(iter=sample_number, burn=burn)    
    return ratio.trace()

### @export 'beta-binomial-model'
alpha = mc.Uninformative('alpha', value=4.)
beta = mc.Uninformative('beta', value=1000.)
pi_mean = mc.Lambda('pi_mean', lambda alpha=alpha, beta=beta: alpha/(alpha+beta))
pi = mc.Beta('pi', alpha, beta, value=pl.maximum(1.e-12, pl.minimum(1-1.e-12, r)))

@mc.potential
def obs(pi=pi):
    return mc.binomial_like(r*n, n, pi)

@mc.deterministic
def pred(alpha=alpha, beta=beta):
    return mc.rbinomial(n_pred, mc.rbeta(alpha, beta)) / float(n_pred)

### @export 'beta-binomial-fit'
mcmc = mc.MCMC([alpha, beta, pi_mean, pi, obs, pred])
mcmc.use_step_method(mc.AdaptiveMetropolis, [alpha, beta])
mcmc.use_step_method(mc.AdaptiveMetropolis, pi)
mcmc.sample(iter*10, burn*10, thin*10)

### @export 'beta-binomial-store'
results['Beta binomial'] = dict(pi=pi_mean.stats(), pred=pred.stats())


### @export 'negative-binomial-model'
pi = mc.Uniform('pi', lower=0, upper=1, value=.5)
#delta = mc.Uninformative('delta', value=100.)
mu_log_10_delta = 1.
log_10_delta = mc.Normal('log_10_delta', mu=mu_log_10_delta, tau=.25**-2)
delta = mc.Lambda('delta', lambda x=log_10_delta: .5 + 10**x)

@mc.potential

import pymc
import radon_varying_intercept
from pylab import hist, show
from pymc import Matplot

M = pymc.MCMC(radon_varying_intercept)
M.sample(iter=50e3, burn=10e3, thin=5)

fit = M.stats()
for k in fit.keys():
     print(k,fit[k]['mean'])

Matplot.plot(M)

def estimate_alpha_beta(data, n_iter=1000, plot=False):
#     data = data.dropna()
    alpha_var = pm.Exponential('alpha', .5)
    beta_var = pm.Exponential('beta', .5)

    observations = pm.Beta('observations', alpha_var, beta_var, value=data,
                           observed=True)

    model = pm.Model([alpha_var, beta_var])
    mcmc = pm.MCMC(model)
    mcmc.sample(n_iter)

    if plot:
        from pymc.Matplot import plot
        plot(mcmc)
        sns.despine()

    alphas = mcmc.trace('alpha')[:]
    betas = mcmc.trace('beta')[:]

    mean_alpha = alphas.mean()
    mean_beta = betas.mean()
    return mean_alpha, mean_beta

#------------------------------------------------------------
# Go through models; compute and plot likelihoods
models = [M0, M1, M2]
linestyles = [':', '--', '-']
labels = ['no outlier correction\n(dotted fit)',
          'mixture model\n(dashed fit)',
          'outlier rejection\n(solid fit)']

x = np.linspace(0, 350, 10)

bins = [(np.linspace(140, 300, 51), np.linspace(0.6, 1.6, 51)),
        (np.linspace(-40, 120, 51), np.linspace(1.8, 2.8, 51)),
        (np.linspace(-40, 120, 51), np.linspace(1.8, 2.8, 51))]

for i, M in enumerate(models):
    S = pymc.MCMC(M)
    S.sample(iter=25000, burn=5000)
    trace = S.trace('beta_M%i' % i)

    H2D, bins1, bins2 = np.histogram2d(trace[:, 0], trace[:, 1], bins=50)
    w = np.where(H2D == H2D.max())

    # choose the maximum posterior slope and intercept
    slope_best = bins1[w[0][0]]
    intercept_best = bins2[w[1][0]]

    # plot the best-fit line
    ax1.plot(x, intercept_best + slope_best * x, linestyles[i], c='k')

    # For the model which identifies bad points,
    # plot circles around points identified as outliers.
    if i == 2:

Date:
    9/22/2013
"""


import pymc
import simple_mixture_model
import sys
import os
sys.path.insert(0, os.path.abspath('..'))
from pysmc import *
import matplotlib.pyplot as plt


mcmc_sampler = pymc.MCMC(simple_mixture_model, verbose=1)
smc_sampler = SMC(mcmc_sampler=mcmc_sampler, num_particles=1000,
                  num_mcmc=10,
                  verbose=3)
smc_sampler.initialize(0.001)
smc_sampler.move_to(1.)
w = smc_sampler.weights
r = smc_sampler.get_particles_of('mixture')
#print w
plt.hist(r, bins=100, weights=w, normed=True)
plt.show()

import warnings
warnings.filterwarnings("ignore")
import diffusion_inverse_model
import matplotlib.pyplot as plt
import sys
import os
sys.path.insert(0, os.path.abspath('..'))
import pysmc as ps
import pymc as pm
import mpi4py.MPI as mpi


if __name__ == '__main__':
    model = diffusion_inverse_model.make_model()
    mcmc = pm.MCMC(model)
    mcmc.use_step_method(ps.GaussianMixtureStep, model['location'])
    mcmc.use_step_method(ps.LognormalRandomWalk, model['alpha'],
                         proposal_sd=1.)
    mcmc.use_step_method(ps.LognormalRandomWalk, model['beta'],
                         proposal_sd=1.)
    mcmc.use_step_method(ps.LognormalRandomWalk, model['tau'],
                         proposal_sd=1.)
    try:
        smc_sampler = ps.SMC(mcmc, num_particles=64, num_mcmc=1,
                             verbose=3, mpi=mpi,
                             db_filename='diffusion_inverse.pickle',
                             update_db=True,
                             gamma_is_an_exponent=True)
        smc_sampler.initialize(0.)
        smc_sampler.move_to(1.)
    except Exception as e:

def compute_MCMC_models(Niter=10000, burn=1000, rseed=0):
    pymc.numpy.random.seed(rseed)

    S1 = pymc.MCMC(model1)
    S1.sample(iter=Niter, burn=burn)
    trace1 = np.vstack([S1.trace('M1_mu')[:],
                        S1.trace('M1_sigma')[:]])
    logp1 = get_logp(S1, model1)

    S2 = pymc.MCMC(model2)
    S2.sample(iter=Niter, burn=burn)
    trace2 = np.vstack([S2.trace('M2_mu1')[:],
                        S2.trace('M2_mu2')[:],
                        S2.trace('M2_sigma1')[:],
                        S2.trace('M2_sigma2')[:],
                        S2.trace('M2_ratio')[:]])
    logp2 = get_logp(S2, model2)

    return trace1, logp1, trace2, logp2

import pymc
import simple
from numpy import mean

model=pymc.MCMC(simple)
model.sample(iter=1000, burn=500, thin=2)

print 'mu',mean(model.trace('mu')[:])
print 'tau',mean(model.trace('tau')[:])