!date
Tue Dec 16 07:55:23 PST 2014
from pymc import *
from scipy.stats.stats import pearsonr
import numpy as np
# set random seed for reproducibility
np.random.seed(12345)
x = np.array([525, 300, 450, 300, 400, 500, 550, 125, 300, 400, 500, 550])
y = np.array([250, 225, 275, 350, 325, 375, 450, 400, 500, 550, 600, 525])
x_mean= x.mean()
y_mean= y.mean()
N = len(x)
pearsonr(x, y)
(0.17466806641204655, 0.58716519218300223)
def model():
data = np.array([x, y])
print(data)
mu1 = Normal('mu1', 0, 0.001)
mu2 = Normal('mu2', 0, 0.001)
lambda1 = Gamma('lambda1', 0.001, 0.001)
lambda2 = Gamma('lambda2', 0.001, 0.001)
rho = Uniform('r', -1, 1)
@pymc.deterministic
def mean(mu1=mu1, mu2=mu2):
return np.array([mu1, mu2])
@pymc.deterministic
def precision(lambda1=lambda1, lambda2=lambda2, rho=rho):
sigma1 = 1 / sqrt(lambda1)
sigma2 = 1 / sqrt(lambda2)
ss1 = sigma1 * sigma2
ss2 = sigma2 * sigma2
rss = rho * sigma1 * sigma2
return [[ss1, rss], [rss, ss2]]
#@pymc.observed
#def obs(value=data, mean=mean, precision=precision):
# return sum(mv_normal_like(v, m, T) for v, m, T in zip(data, mean, precision))
xy = MvNormal('xy', mu=mean, tau=precision, value=data.T, observed=True)
M = pymc.MCMC(locals())
M.sample(10000, 5000)
return M
M = model()
[[525 300 450 300 400 500 550 125 300 400 500 550] [250 225 275 350 325 375 450 400 500 550 600 525]] [-----------------100%-----------------] 10000 of 10000 complete in 38.4 sec
/homes/abie/anaconda/lib/python2.7/site-packages/pymc/NumpyDeterministics.py:187: RuntimeWarning: divide by zero encountered in power sqrt_jacobians = {'x': lambda x: .5 * x ** -.5} /homes/abie/anaconda/lib/python2.7/site-packages/pymc/CommonDeterministics.py:671: RuntimeWarning: divide by zero encountered in double_scalars return op_function_base(a, b) /homes/abie/anaconda/lib/python2.7/site-packages/pymc/CommonDeterministics.py:828: RuntimeWarning: divide by zero encountered in divide truediv_jacobians = {'a': lambda a, b: ones(shape(a)) / b, /homes/abie/anaconda/lib/python2.7/site-packages/pymc/CommonDeterministics.py:829: RuntimeWarning: divide by zero encountered in double_scalars 'b': lambda a, b: - a / b ** 2}