%reset -f
from pymc import *
from scipy.stats.stats import pearsonr
import numpy as np
%matplotlib inline
print(pymc.__version__)
2.3.3
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])
N = len(x)
data = np.array([x, y])
mean = data.mean(1)
std = data.std(1)
print(data)
print(mean)
print(std)
[[525 300 450 300 400 500 550 125 300 400 500 550] [250 225 275 350 325 375 450 400 500 550 600 525]] [ 408.33333333 402.08333333] [ 125.13881181 118.34727031]
pearsonr(x, y)
(0.1746680664120466, 0.58716519218300223)
mu1 = Normal('mu1', 0, 0.001, value=mean[0], observed=True)
mu2 = Normal('mu2', 0, 0.001, value=mean[1], observed=True)
lambda1 = Gamma('lambda1', 0.001, 0.001)
lambda2 = Gamma('lambda2', 0.001, 0.001)
rho = Uniform('rho', -1, 1, value=0)
@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 * sigma1
ss2 = sigma2 * sigma2
rss = rho * sigma1 * sigma2
#return np.power(np.mat([[ss1, rss], [rss, ss2]]), -1)
return np.mat([[ss1, rss], [rss, ss2]])
xy = MvNormal('xy', mu=mean, tau=precision, value=data.T, observed=True)
M = pymc.MCMC(locals())
M.sample(100, 50)
M.db.close()
[-----------------100%-----------------] 100 of 100 complete in 0.3 sec
/Library/Python/2.7/site-packages/pymc/NumpyDeterministics.py:188: RuntimeWarning: divide by zero encountered in power sqrt_jacobians = {'x': lambda x: .5 * x ** -.5} /Library/Python/2.7/site-packages/pymc/CommonDeterministics.py:682: RuntimeWarning: divide by zero encountered in double_scalars return op_function_base(a, b) /Library/Python/2.7/site-packages/pymc/CommonDeterministics.py:843: RuntimeWarning: divide by zero encountered in divide truediv_jacobians = {'a': lambda a, b: ones(shape(a)) / b, /Library/Python/2.7/site-packages/pymc/CommonDeterministics.py:844: RuntimeWarning: divide by zero encountered in double_scalars 'b': lambda a, b: - a / b ** 2} /Library/Python/2.7/site-packages/pymc/CommonDeterministics.py:682: RuntimeWarning: invalid value encountered in multiply return op_function_base(a, b)
pymc.Matplot.plot(M)
Plotting mean_0 Plotting mean_1 Plotting lambda1 Plotting rho Plotting
/usr/local/lib/python2.7/site-packages/numpy/core/fromnumeric.py:2507: VisibleDeprecationWarning: `rank` is deprecated; use the `ndim` attribute or function instead. To find the rank of a matrix see `numpy.linalg.matrix_rank`. VisibleDeprecationWarning) /usr/local/lib/python2.7/site-packages/matplotlib/axes/_axes.py:1743: RuntimeWarning: divide by zero encountered in true_divide c /= np.sqrt(np.dot(x, x) * np.dot(y, y))
rho.get_value()
array(-0.9323782802284845)