{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 1200x1200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 1200x1200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 1200x1200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 1200x1200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "for alpha in [[1, 3, 4], [1, 1, 1], [10, 0.2, 0.2], [0.1, 0.1, 0.1]]:\n",
    "    theta = stats.dirichlet(alpha).rvs(500)\n",
    "\n",
    "    fig = plt.figure(figsize=(8, 8), dpi=150)\n",
    "    ax = plt.gca(projection='3d')\n",
    "    plt.title(r'$\\alpha$ = {}'.format(alpha))\n",
    "    ax.scatter(theta[:, 0], theta[:, 1], theta[:, 2])\n",
    "    ax.view_init(azim=30)\n",
    "    ax.set_xlabel(r'$\\theta_1$')\n",
    "    ax.set_ylabel(r'$\\theta_2$')\n",
    "    ax.set_zlabel(r'$\\theta_3$')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "from sklearn.datasets import load_iris\n",
    "import seaborn as sns\n",
    "from scipy import stats\n",
    "import tensorflow as tf\n",
    "import edward as ed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(load_iris()['data'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<seaborn.axisgrid.PairGrid at 0x7f92097808d0>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 93140 characters]\n",
      "text/plain": [
       "<Figure size 762.375x720 with 20 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df['target'] = load_iris()['target']\n",
    "sns.pairplot(df, hue='target', vars=[0, 1, 2, 3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 4.,  1.,  6.,  5., 12.,  8.,  4.,  5.,  2.,  3.]),\n",
       " array([4.3 , 4.45, 4.6 , 4.75, 4.9 , 5.05, 5.2 , 5.35, 5.5 , 5.65, 5.8 ]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 12684 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(12, 6))\n",
    "b = 8\n",
    "a = 4\n",
    "\n",
    "plt.title('MOAS')\n",
    "plt.xlabel(r'$X$')\n",
    "plt.ylabel('bin count')\n",
    "plt.annotate('Cummulative distribution', (6.8, 17), (7, 20), arrowprops=dict(facecolor='black', shrink=0.05))\n",
    "plt.hist(df[0], color=sns.color_palette()[3], rwidth=a)\n",
    "plt.hist(df[0][df['target'] == 2], color=sns.color_palette()[2], rwidth=a)\n",
    "plt.hist(df[0][df['target'] == 1], color=sns.color_palette()[1], rwidth=a)\n",
    "plt.hist(df[0][df['target'] == 0], color=sns.color_palette()[0], rwidth=a)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f92086cec88>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 52304 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(13)\n",
    "x = np.linspace(-2, 6, 500)\n",
    "pi = stats.dirichlet([1, 1, 1]).rvs()[0]\n",
    "mu = stats.norm(1, 1).rvs(3)\n",
    "\n",
    "y = stats.norm(mu, np.ones(3)).pdf(x[:, None])\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.title('Gaussian Mixture model')\n",
    "plt.ylabel(r'$P(x)$')\n",
    "plt.xlabel(r'$X$')\n",
    "\n",
    "plt.plot(x, y[:, 1] * pi[1], ls='--')\n",
    "plt.plot(x, y[:, 2] * pi[2], ls='--')\n",
    "plt.plot(x, y[:, 0] * pi[0], ls='--')\n",
    "\n",
    "for i in range(3):\n",
    "    xi = x[np.argmax(y[:, i])]\n",
    "    yi = (y[:, i] * pi[i]).max()\n",
    "    plt.vlines(xi, 0, yi)\n",
    "    plt.text(xi + 0.05, yi + 0.01, r'$\\pi_{}$'.format(i + 1))\n",
    "plt.plot(x, (y * pi).sum(1))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(load_iris()['data'])\n",
    "y = df.values\n",
    "y = (y - y.mean(axis=0)) / y.std(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 3 # number of clusters\n",
    "\n",
    "d = df.shape[1]\n",
    "n = df.shape[0]\n",
    "\n",
    "pi = ed.models.Dirichlet(tf.ones(k))\n",
    "mu = ed.models.Normal(tf.zeros(d), tf.ones(d), sample_shape=k)  # shape (3, 4) 3 gaussians, 4 variates\n",
    "sigmasq = ed.models.InverseGamma(tf.ones(d), tf.ones(d), sample_shape=k)\n",
    "x = ed.models.ParamMixture(pi, {'loc': mu, 'scale_diag': tf.sqrt(sigmasq)},\n",
    "                 ed.models.MultivariateNormalDiag,\n",
    "                 sample_shape=n)\n",
    "z = x.cat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = 500  # number of samples\n",
    "\n",
    "qpi = ed.models.Empirical(tf.get_variable('qpi', shape=[t, k], initializer=tf.constant_initializer(1 / k)))\n",
    "qmu = ed.models.Empirical(tf.get_variable('qmu', shape=[t, k, d], initializer=tf.zeros_initializer()))\n",
    "qsigmasq = ed.models.Empirical(tf.get_variable('qsigmasq', shape=[t, k, d], initializer=tf.ones_initializer()))\n",
    "qz = ed.models.Empirical(tf.get_variable('qz', shape=[t, n], initializer=tf.zeros_initializer(), dtype=tf.int32))     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.6/site-packages/edward/util/random_variables.py:52: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  not np.issubdtype(value.dtype, np.float) and \\\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "500/500 [100%] ██████████████████████████████ Elapsed: 2s | Acceptance Rate: 1.000\n"
     ]
    }
   ],
   "source": [
    "inference = ed.Gibbs({pi: qpi, mu: qmu, sigmasq: qsigmasq, z: qz}, \n",
    "                    data={x: y})\n",
    "inference.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_s = qmu.sample(500).eval()\n",
    "sigmasq_s = qsigmasq.sample(500).eval()\n",
    "pi_s = qpi.sample(500).eval()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 57488 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(12, 6))\n",
    "plt.title(r'Distribution of the likelihood of $\\mu | X$ in all dimensions')\n",
    "plt.ylabel(r'$P(\\mu|X)$')\n",
    "plt.xlabel(r'$\\mu$')\n",
    "for i in range(3):\n",
    "    for j in range(4):\n",
    "        sns.distplot(mu_s[:, i, j], hist=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,\n",
       "       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,\n",
       "       2, 2, 2, 2, 2, 2, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n",
       "       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1,\n",
       "       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0,\n",
       "       0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0,\n",
       "       0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1])"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.vstack([pi_s.mean(0)[i] * \\\n",
    "           stats.multivariate_normal(mu_s.mean(0)[i], np.sqrt(sigmasq_s.mean(0))[i]).pdf(y) \\\n",
    "           for i in range(3)]).argmax(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.94"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.vstack([pi_s.mean(0)[i] * \\\n",
    "           stats.multivariate_normal(mu_s.mean(0)[i], np.sqrt(sigmasq_s.mean(0))[i]).pdf(y) \\\n",
    "          for i in range(3)])[::-1].argmax(axis=0).sum() / df.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<BarContainer object of 3 artists>"
      ]
     },
     "execution_count": 114,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 6536 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar(np.arange(pi_s.shape[1]), pi_s.mean(0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 188,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 23064 characters]\n",
      "text/plain": [
       "<Figure size 1008x288 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 22952 characters]\n",
      "text/plain": [
       "<Figure size 1008x288 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23596 characters]\n",
      "text/plain": [
       "<Figure size 1008x288 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 27036 characters]\n",
      "text/plain": [
       "<Figure size 1008x288 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 30748 characters]\n",
      "text/plain": [
       "<Figure size 1008x288 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def dirichlet_process(h_0, alpha):\n",
    "    \"\"\"\n",
    "    Truncated dirichlet process.\n",
    "    :param h_0: (scipy distribution)\n",
    "    :param alpha: (flt)\n",
    "    :param n: (int) Truncate value.\n",
    "    \"\"\"\n",
    "    n = max(int(5 * alpha + 2), 500)  # truncate the values. \n",
    "    pi = stats.beta(1, alpha).rvs(size=n)\n",
    "    pi[1:] = pi[1:] * (1 - pi[:-1]).cumprod()  # stick breaking process\n",
    "    theta = h_0(size=n)  # samples from original distribution\n",
    "    return pi, theta\n",
    "        \n",
    "def plot_normal_dp_approximation(alpha, n=3):\n",
    "    pi, theta = dirichlet_process(stats.norm.rvs, alpha)\n",
    "    x = np.linspace(-4, 4, 100)\n",
    "    \n",
    "    plt.figure(figsize=(14, 4))\n",
    "    plt.suptitle(r'{} samples from DP($\\alpha$). $\\alpha$ = {}'.format(n, alpha))\n",
    "    plt.ylabel(r'$\\pi$')\n",
    "    plt.xlabel(r'$\\theta$')\n",
    "    pltcount = int('1' + str(n) + '0')\n",
    "    \n",
    "    for i in range(n):\n",
    "        pltcount += 1\n",
    "        plt.subplot(pltcount)\n",
    "        pi, theta = dirichlet_process(stats.norm.rvs, alpha)\n",
    "        pi = pi * (stats.norm.pdf(0) / pi.max())\n",
    "        plt.vlines(theta, 0, pi, alpha=0.5)\n",
    "        plt.ylim(0, 1)\n",
    "        plt.plot(x, stats.norm.pdf(x))\n",
    "\n",
    "np.random.seed(65)\n",
    "for alpha in [.1, 1, 10, 100, 1000]:\n",
    "    plot_normal_dp_approximation(alpha)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.LineCollection at 0x7f91bfd9bfd0>"
      ]
     },
     "execution_count": 195,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 47556 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(95)\n",
    "\n",
    "acceptance_p = 0.98\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.title(f'Number samples need to have a total probability equal to {acceptance_p}%')\n",
    "plt.ylabel('$\\sum_{i=1}^{\\infty}{\\pi_i}$')\n",
    "plt.xlabel('Number of clusters $k$')\n",
    "def plot_summation(alpha):\n",
    "    pi, _ = dirichlet_process(stats.norm.rvs, alpha)\n",
    "    p_total = np.cumsum(pi)\n",
    "    i = np.argmin(np.abs(p_total - acceptance_p))\n",
    "    plt.plot(np.arange(pi[:i + 1].shape[0]), p_total[:i + 1])\n",
    "    return i\n",
    "\n",
    "k = 0\n",
    "for alpha in [1, 2, 5, 7, 10]:\n",
    "    k_i = plot_summation(alpha)\n",
    "    k = max(k, k_i)\n",
    "    plt.vlines(k_i, 0, acceptance_p, linestyles='--')\n",
    "    plt.text(k_i + 0.2, 0.2, f'k = {k_i}', rotation=45)\n",
    "    \n",
    "plt.ylim(0, 1)\n",
    "plt.xlim(0, k + 2)\n",
    "plt.hlines(acceptance_p, 0, k, linestyles='--')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
      "  warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f20ba3e1828>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 11660 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(12, 6))\n",
    "plt.title('Histogram of the 3d column of the (standardized) Iris dataset.')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('count')\n",
    "sns.distplot(y[:, 2], bins=20, kde=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymc3 as pm\n",
    "from theano import tensor as tt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def stick_breaking(beta):\n",
    "    portion_remaining = tt.concatenate([[1], tt.extra_ops.cumprod(1 - beta)[:-1]])\n",
    "    return beta * portion_remaining"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f20b78504e0>]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 13100 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0, 1, 100)\n",
    "plt.plot(x, stats.gamma(1, 0, 0.5).pdf(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f20b78300f0>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 12884 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x, stats.halfnorm(0, 1).pdf(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0295230814623677"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.gamma(1, 0, 1).rvs(1000).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 35\n",
    "m = pm.Model()\n",
    "with m:\n",
    "    alpha = pm.Gamma('alpha', 1, 1)\n",
    "    beta = pm.Beta('beta', 1, alpha, shape=k)\n",
    "    pi = pm.Deterministic('pi', stick_breaking(beta))\n",
    "\n",
    "    sd = pm.HalfNormal('sd', 1, shape=k)    \n",
    "    lambda_1 = pm.Uniform('lambda1', 0, 10, shape=k)\n",
    "    mu_0 = pm.Normal('mu_0', 0, 1)\n",
    "    \n",
    "    mu = pm.Normal('mu', mu=mu_0, sd=lambda_1 * sd, shape=k)\n",
    "    obs = pm.NormalMixture('obs', w=pi, mu=mu, sd=lambda_ * sd, observed=y[:, 2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 35\n",
    "m = pm.Model()\n",
    "with m:\n",
    "    alpha = pm.Gamma('alpha', 1, 1)\n",
    "    beta = pm.Beta('beta', 1, alpha, shape=k)\n",
    "    pi = pm.Deterministic('pi', stick_breaking(beta))\n",
    "\n",
    "    tau = pm.Gamma('tau', 1, 1, shape=k)    \n",
    "    lambda_ = pm.Uniform('lambda', 0, 1, shape=k)\n",
    "    mu0 = pm.Uniform('mu0', -3, 3, shape=k)\n",
    "    \n",
    "    mu = pm.Normal('mu', mu=mu0, tau=1, shape=k)\n",
    "    obs = pm.NormalMixture('obs', w=pi, mu=mu, tau=lambda_ * tau, observed=y[:, 2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha_log__ -1.0\n",
      "beta_logodds__ -48.52030263919613\n",
      "tau_log__ -35.0\n",
      "lambda1_interval__ -48.52030263919613\n",
      "lambda2 -56.42299998176168\n",
      "mu -16.127760854365846\n",
      "obs -256.61897509450506\n"
     ]
    }
   ],
   "source": [
    "for RV in m.basic_RVs:\n",
    "    print(RV.name, RV.logp(m.test_point))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 245,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using advi...\n",
      "Average Loss = 265.37:   8%|▊         | 15305/200000 [00:13<02:46, 1106.63it/s]\n",
      "Convergence archived at 15400\n",
      "Interrupted at 15,399 [7%]: Average Loss = 275.51\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [mu, mu0_interval__, lambda_interval__, tau_log__, beta_logodds__, alpha_log__]\n",
      "100%|██████████| 1000/1000 [02:09<00:00,  7.71it/s]\n",
      "There were 22 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "The acceptance probability does not match the target. It is 0.8832672993179058, but should be close to 0.8. Try to increase the number of tuning steps.\n",
      "There were 60 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "The acceptance probability does not match the target. It is 0.9412707023823474, but should be close to 0.8. Try to increase the number of tuning steps.\n",
      "There were 118 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "The acceptance probability does not match the target. It is 0.8960563359041548, but should be close to 0.8. Try to increase the number of tuning steps.\n",
      "There were 41 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "The gelman-rubin statistic is larger than 1.4 for some parameters. The sampler did not converge.\n",
      "The estimated number of effective samples is smaller than 200 for some parameters.\n"
     ]
    }
   ],
   "source": [
    "with m: \n",
    "    step = None\n",
    "    trace = pm.sample(500, tune=500, init='advi', random_seed=35171, step=step)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[<matplotlib.axes._subplots.AxesSubplot object at 0x7f2046b33048>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7f2046adaba8>],\n",
       "       [<matplotlib.axes._subplots.AxesSubplot object at 0x7f2046a84e80>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7f2046aa6e48>],\n",
       "       [<matplotlib.axes._subplots.AxesSubplot object at 0x7f2046a52198>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7f20469f7438>],\n",
       "       [<matplotlib.axes._subplots.AxesSubplot object at 0x7f2046a1f710>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7f20469ca9e8>],\n",
       "       [<matplotlib.axes._subplots.AxesSubplot object at 0x7f20469f2cc0>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7f204699cf98>],\n",
       "       [<matplotlib.axes._subplots.AxesSubplot object at 0x7f204694f2b0>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7f20468f8588>],\n",
       "       [<matplotlib.axes._subplots.AxesSubplot object at 0x7f2046921860>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7f20468ccb38>]],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 120,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 595336 characters]\n",
      "text/plain": [
       "<Figure size 864x1008 with 14 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.traceplot(trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 12560 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(12, 6))\n",
    "plt.title(r\"Distribution of $\\pi'$\")\n",
    "plt.bar(np.arange(k) + 1 - 0.5, trace['pi'].mean(axis=0), width=1., lw=0);\n",
    "plt.xlim(0.5, k);\n",
    "plt.xlabel('Component');\n",
    "plt.ylabel('Posterior expected mixture weight');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 189,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 30188 characters]\n",
      "text/plain": [
       "<Figure size 1008x288 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_normal_dp_approximation(trace['alpha'].mean(), 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 210,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.64885315, -1.30709587,  0.02895939, -0.04105907, -0.08588051,\n",
       "       -0.08249138, -0.13069356, -0.06956456, -0.11705567, -0.28732056,\n",
       "       -0.14913065, -0.09560382, -0.17197101, -0.18029005, -0.13289487,\n",
       "       -0.0966986 , -0.14822977, -0.15511884, -0.01362059, -0.08300809,\n",
       "       -0.12929974, -0.16073243, -0.16999157, -0.19414815, -0.12453189,\n",
       "       -0.15904975, -0.10099628, -0.0359372 , -0.14126203, -0.18965758,\n",
       "       -0.19326695, -0.11220963, -0.10462724, -0.2159177 , -0.19433511])"
      ]
     },
     "execution_count": 210,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 243,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
      "  warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f205b91be10>"
      ]
     },
     "execution_count": 243,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 26680 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(12, 6))\n",
    "plt.title('Density estimation with Dirichlet Processes')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('p(x)')\n",
    "x = np.linspace(-3,  3, 100)\n",
    "pi_mean = trace['pi'].mean(0)\n",
    "mu_mean = trace['mu'].mean(0) \n",
    "mu0_mean = trace['mu0'].mean(0)\n",
    "lambda_mean = trace['lambda1'].mean(0)\n",
    "tau_mean = trace['tau'].mean(0)\n",
    "sd_mean = np.sqrt((1 / tau_mean * lambda_mean))\n",
    "tau_sd = trace['mu'].std(0)\n",
    "lambda_sd = trace['lambda1'].std(0)\n",
    "sd_sd = np.sqrt((1 / tau_sd * lambda_sd))\n",
    "\n",
    "p = stats.norm.interval(0.05, trace['mu'].mean(0), trace['mu'].std(0))\n",
    "psd = stats.norm.interval(0.05, sd_mean, sd_std)\n",
    "\n",
    "a = 2\n",
    "plt.fill_between(x, stats.norm.pdf(x, mu_mean[0], psd[0][0]) * pi_mean[0], stats.norm.pdf(x, mu_mean[0], psd[1][0]) * pi_mean[0])\n",
    "plt.fill_between(x, stats.norm.pdf(x, mu_mean[1], psd[0][1]) * pi_mean[1], stats.norm.pdf(x, mu_mean[1], psd[1][1]) * pi_mean[1])\n",
    "sns.distplot(y[:, 2], bins=20, kde=None, norm_hist=True, rug=True, )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
      "  warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 25460 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(12, 6))\n",
    "plt.title('Density estimation with Dirichlet Processes')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('p(x)')\n",
    "x = np.linspace(-3,  3, 100)\n",
    "pi_mean = trace['pi'].mean(0)\n",
    "mu_mean = trace['mu'].mean(0) \n",
    "mu0_mean = trace['mu0'].mean(0)\n",
    "lambda_1_mean = trace['lambda1'].mean(0)\n",
    "tau_mean = trace['tau'].mean(0)\n",
    "sd_mean = np.sqrt((1 / tau_mean * lambda_1_mean))\n",
    "\n",
    "plt.plot(x, stats.norm.pdf(np.repeat(x[:, None], k, axis=1), \n",
    "                           loc=mu_mean, \n",
    "                           scale=sd_mean) * pi_mean,\n",
    "        ls='--')\n",
    "sns.distplot(y[:, 2], bins=20, kde=None, norm_hist=True, rug=True, )\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
