{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.3.3\n"
     ]
    }
   ],
   "source": [
    "%reset -f\n",
    "from pymc import *\n",
    "from scipy.stats.stats import pearsonr\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "print(pymc.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[525 300 450 300 400 500 550 125 300 400 500 550]\n",
      " [250 225 275 350 325 375 450 400 500 550 600 525]]\n",
      "[ 408.33333333  402.08333333]\n",
      "[ 125.13881181  118.34727031]\n"
     ]
    }
   ],
   "source": [
    "x = np.array([525, 300, 450, 300, 400, 500, 550, 125, 300, 400, 500, 550])\n",
    "y = np.array([250, 225, 275, 350, 325, 375, 450, 400, 500, 550, 600, 525])\n",
    "N = len(x)\n",
    "data = np.array([x, y])\n",
    "mean = data.mean(1)\n",
    "std = data.std(1)\n",
    "print(data)\n",
    "print(mean)\n",
    "print(std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.1746680664120466, 0.58716519218300223)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pearsonr(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r",
      " [-----------------100%-----------------] 100 of 100 complete in 0.3 sec"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Python/2.7/site-packages/pymc/NumpyDeterministics.py:188: RuntimeWarning: divide by zero encountered in power\n",
      "  sqrt_jacobians = {'x': lambda x: .5 * x ** -.5}\n",
      "/Library/Python/2.7/site-packages/pymc/CommonDeterministics.py:682: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  return op_function_base(a, b)\n",
      "/Library/Python/2.7/site-packages/pymc/CommonDeterministics.py:843: RuntimeWarning: divide by zero encountered in divide\n",
      "  truediv_jacobians = {'a': lambda a, b: ones(shape(a)) / b,\n",
      "/Library/Python/2.7/site-packages/pymc/CommonDeterministics.py:844: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  'b': lambda a, b: - a / b ** 2}\n",
      "/Library/Python/2.7/site-packages/pymc/CommonDeterministics.py:682: RuntimeWarning: invalid value encountered in multiply\n",
      "  return op_function_base(a, b)\n"
     ]
    }
   ],
   "source": [
    "mu1 = Normal('mu1', 0, 0.001, value=mean[0], observed=True)\n",
    "mu2 = Normal('mu2', 0, 0.001, value=mean[1], observed=True)\n",
    "lambda1 = Gamma('lambda1', 0.001, 0.001)\n",
    "lambda2 = Gamma('lambda2', 0.001, 0.001)\n",
    "rho = Uniform('rho', -1, 1, value=0)\n",
    "\n",
    "@pymc.deterministic\n",
    "def mean(mu1=mu1, mu2=mu2):\n",
    "    return np.array([mu1, mu2])\n",
    "\n",
    "@pymc.deterministic\n",
    "def precision(lambda1=lambda1, lambda2=lambda2, rho=rho):\n",
    "    sigma1 = 1 / sqrt(lambda1)\n",
    "    sigma2 = 1 / sqrt(lambda2)\n",
    "    ss1 = sigma1 * sigma1\n",
    "    ss2 = sigma2 * sigma2\n",
    "    rss = rho * sigma1 * sigma2\n",
    "    #return np.power(np.mat([[ss1, rss], [rss, ss2]]), -1)\n",
    "    return np.mat([[ss1, rss], [rss, ss2]])\n",
    "\n",
    "xy = MvNormal('xy', mu=mean, tau=precision, value=data.T, observed=True)\n",
    "\n",
    "M = pymc.MCMC(locals())\n",
    "M.sample(100, 50)\n",
    "M.db.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Plotting mean_0\n",
      "Plotting mean_1\n",
      "Plotting lambda1\n",
      "Plotting rho\n",
      "Plotting"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/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`.\n",
      "  VisibleDeprecationWarning)\n",
      "/usr/local/lib/python2.7/site-packages/matplotlib/axes/_axes.py:1743: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  c /= np.sqrt(np.dot(x, x) * np.dot(y, y))\n"
     ]
    }
   ],
   "source": [
    "pymc.Matplot.plot(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(-0.9323782802284845)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rho.get_value()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}