#!/usr/bin/env python
# coding: utf-8

# In[1]:


get_ipython().run_line_magic('reset', '-f')
from pymc import *
from scipy.stats.stats import pearsonr
import numpy as np
get_ipython().run_line_magic('matplotlib', 'inline')
print(pymc.__version__)


# In[2]:


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)


# In[3]:


pearsonr(x, y)


# In[4]:


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()


# In[ ]:


pymc.Matplot.plot(M)


# In[9]:


rho.get_value()


# In[ ]:





# In[ ]: