MCMC multitau.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import pymc3 as pm\n",
    "import theano.tensor as tt\n",
    "\n",
    "\n",
    "SMALL_SIZE = 16\n",
    "MEDIUM_SIZE = 18\n",
    "BIGGER_SIZE = 20\n",
    "\n",
    "plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "\n",
    "SEED = 35010732 # from random.org\n",
    "np.random.seed(SEED)\n",
    "\n",
    "print(plt.style.available)\n",
    "plt.style.use('seaborn-white')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The conditional probability for a OU process $p(x,t|x_{0},0)$ is\n",
    "\n",
    "$$p(x,t|x_{0},0)=\\frac{1}{\\sqrt{2\\pi A(1-B^{2})}}\\exp \\left(-\\frac{(x-Bx_{0})^{2}}{2A(1-B^{2})}\\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Ornstein_Uhlenbeck(pm.Continuous):\n",
    "    \"\"\"\n",
    "    Ornstein-Uhlenbeck Process\n",
    "    Parameters\n",
    "    ----------\n",
    "    B : tensor\n",
    "        B > 0, B = exp(-(D/A)*delta_t)\n",
    "    A : tensor\n",
    "        A > 0, amplitude of fluctuation <x**2>=A\n",
    "    delta_t: scalar\n",
    "        delta_t > 0, time step\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, A=None, B=None,\n",
    "                 *args, **kwargs):\n",
    "        super(Ornstein_Uhlenbeck, self).__init__(*args, **kwargs)\n",
    "        self.A = A\n",
    "        self.B = B\n",
    "        self.mean = 0.\n",
    "\n",
    "    def logp(self, x):\n",
    "        A = self.A\n",
    "        B = self.B\n",
    "\n",
    "        x_im1 = x[:-1]\n",
    "        x_i = x[1:]\n",
    "\n",
    "        ou_like = pm.Normal.dist(mu=x_im1*B, tau=1.0/A/(1-B**2)).logp(x_i)\n",
    "        return pm.Normal.dist(mu=0.0,tau=1.0/A).logp(x[0]) + tt.sum(ou_like)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.load(\"OUmt_sN05.npy\")\n",
    "data = data[:2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a_bound = 10\n",
    "result_df = pd.DataFrame(columns=['dt','A', 'dA','B','dB','s','ds'])\n",
    "for dataset in data:\n",
    "    delta_t = dataset[0]\n",
    "    ts = dataset[1:]\n",
    "    print(delta_t)\n",
    "    with pm.Model() as model:\n",
    "        B = pm.Beta('B', alpha=1.0,beta=1.0)\n",
    "        A = pm.Uniform('A', lower=0, upper=a_bound)\n",
    "        sigma = pm.Uniform('sigma',lower=0,upper=5)\n",
    "\n",
    "        path = Ornstein_Uhlenbeck('path',A=A, B=B,shape=len(ts))\n",
    "        dataObs = pm.Normal('dataObs',mu=path,sigma=sigma,observed=ts)\n",
    "        trace = pm.sample(2000,cores=4)\n",
    "        \n",
    "    a_mean = trace['A'].mean()\n",
    "    b_mean = trace['B'].mean()\n",
    "    a_std = trace['A'].std()\n",
    "    b_std = trace['B'].std()\n",
    "    sigma_mean = trace['sigma'].mean()\n",
    "    sigma_std = trace['sigma'].std()\n",
    "    \n",
    "    result_df = result_df.append({'dt':delta_t,\n",
    "                      'A':a_mean,\n",
    "                      'dA':a_std,\n",
    "                      'B':b_mean,\n",
    "                      'dB':b_std,\n",
    "                      's':sigma_mean,\n",
    "                      'ds':sigma_std},ignore_index=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tau = -delta_t_list/np.log(result_array.T[2])\n",
    "dtau = delta_t_list*result_array.T[3]/result_array.T[2]/np.log(result_array.T[2])**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(delta_t_list,result_array.T[6],\"o\")\n",
    "plt.xlabel(r'$\\Delta t/\\tau$')\n",
    "plt.ylabel(r'$\\sigma_{GT-model}$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(delta_t_list,result_array.T[0],yerr=result_array.T[1],fmt=\"o\",label=\"A\")\n",
    "plt.errorbar(delta_t_list,tau,dtau,fmt=\"o\",label=r'$\\tau$')\n",
    "plt.legend(loc=\"upper left\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(delta_t_list,result_array.T[4],yerr=result_array.T[5],fmt=\"o\")\n",
    "plt.xlabel(r'$\\Delta t/\\tau$')\n",
    "plt.ylabel(r'$\\sigma_{noise}$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}