!date

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)

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