{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pymc3 as pm\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "from scipy import stats\n",
    "import pandas as pd\n",
    "import theano\n",
    "import theano.tensor as tt\n",
    "from fbprophet import Prophet\n",
    "import seaborn as sns\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(style='white')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7fddaa90bef0>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 68388 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = pd.read_csv('peyton_manning.csv')\n",
    "# Make sure we work with datetime types\n",
    "df['ds'] = pd.to_datetime(df['ds'])\n",
    "# Scale the data\n",
    "df['y_scaled'] = df['y'] / df['y'].max()\n",
    "df['t'] = (df['ds'] - df['ds'].min()) / (df['ds'].max() - df['ds'].min())\n",
    "df.plot(x='ds', y='y', figsize=(16, 6), title=\"Wikipedia page views for 'Peyton Manning'\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Linear trend with changepoints\n",
    "$$ g(t) = (k + a(t)^T \\delta)t + (m + a(t)^T \\gamma) $$\n",
    "\n",
    "$y_j$ is set to $-s_j \\delta_j$\n",
    "\n",
    "$$ \\vec{g} = (k + A \\delta) \\odot t + (m + A \\gamma) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 49660 characters]\n",
      "text/plain": [
       "<Figure size 1152x648 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(25)\n",
    "n_changepoints = 10\n",
    "t = np.arange(1000)\n",
    "s = np.sort(np.random.choice(t, n_changepoints, replace=False))\n",
    "\n",
    "A = (t[:, None] > s) * 1\n",
    "delta = np.random.normal(size=n_changepoints)\n",
    "k = 1\n",
    "m = 5\n",
    "\n",
    "growth = (k + A @ delta) * t\n",
    "gamma = -s * delta\n",
    "offset = m + A @ gamma\n",
    "trend = growth + offset\n",
    "\n",
    "plt.figure(figsize=(16, 3 * 3))\n",
    "n = 310\n",
    "i = 0\n",
    "for t, f in zip(['Linear Trend with Changepoints', 'Growth rate', 'Growth offset'], \n",
    "                [trend, growth, offset]):\n",
    "    i += 1\n",
    "    plt.subplot(n + i)\n",
    "    plt.title(t)\n",
    "    plt.yticks([])\n",
    "    plt.vlines(s, min(f), max(f), lw=0.8, linestyles='--')\n",
    "    plt.plot(f)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def det_dot(a, b):\n",
    "    \"\"\"\n",
    "    The theano dot product and NUTS sampler doen't work with large matrices?\n",
    "    \n",
    "    :param a: (np matrix)\n",
    "    :param b: (theano vector)\n",
    "    \"\"\"\n",
    "    return (a * b[None, :]).sum(axis=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trend_model(m, t, n_changepoints=25, changepoints_prior_scale=0.05, \n",
    "                growth_prior_scale=5, changepoint_range=0.8):\n",
    "    \"\"\"\n",
    "    The piecewise linear trend with changepoint implementation in PyMC3.\n",
    "    :param m: (pm.Model)\n",
    "    :param t: (np.array) MinMax scaled time.\n",
    "    :param n_changepoints: (int) The number of changepoints to model.\n",
    "    :param changepoint_prior_scale: (flt/ None) The scale of the Laplace prior on the delta vector.\n",
    "                                    If None, a hierarchical prior is set.\n",
    "    :param growth_prior_scale: (flt) The standard deviation of the prior on the growth.\n",
    "    :param changepoint_range: (flt) Proportion of history in which trend changepoints will be estimated. \n",
    "    :return g, A, s: (tt.vector, np.array, tt.vector)\n",
    "    \"\"\"\n",
    "    s = np.linspace(0, changepoint_range * np.max(t), n_changepoints + 1)[1:]\n",
    "    \n",
    "    # * 1 casts the boolean to integers\n",
    "    A = (t[:, None] > s) * 1\n",
    "\n",
    "    with m:\n",
    "        # initial growth\n",
    "        k = pm.Normal('k', 0 , growth_prior_scale)\n",
    "        \n",
    "        if changepoints_prior_scale is None:\n",
    "            changepoints_prior_scale = pm.Exponential('tau', 1.5, testval=0.05)\n",
    "        \n",
    "        # rate of change\n",
    "        delta = pm.Laplace('delta', 0, changepoints_prior_scale, shape=n_changepoints)\n",
    "        # offset\n",
    "        m = pm.Normal('m', 0, 5)\n",
    "        gamma = -s * delta\n",
    "        \n",
    "        g = (k + det_dot(A, delta)) * t + (m + det_dot(A, gamma))\n",
    "    return g, A, s\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = pm.Model()\n",
    "\n",
    "with m:\n",
    "    y, A, s = trend_model(m, df['t'])\n",
    "    \n",
    "    sigma = pm.HalfCauchy('sigma', 0.5, testval=1)\n",
    "    pm.Normal('obs', \n",
    "                 mu=y,\n",
    "                 sd=sigma,\n",
    "                 observed=df['y_scaled'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 140532 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def sanity_check(m, df):\n",
    "    \"\"\"\n",
    "    :param m: (pm.Model)\n",
    "    :param df: (pd.DataFrame)\n",
    "    \"\"\"\n",
    "    # Sample from the prior and check of the model is well defined.\n",
    "    y = pm.sample_prior_predictive(model=m, vars=['obs'])['obs']\n",
    "    plt.figure(figsize=(16, 6))\n",
    "    plt.plot(y.mean(0), label='mean prior')\n",
    "    plt.fill_between(np.arange(y.shape[1]), -y.std(0), y.std(0), alpha=0.25, label='standard deviation')\n",
    "    plt.plot(df['y_scaled'], label='true value')\n",
    "    plt.legend()\n",
    "\n",
    "# And run the sanity check\n",
    "sanity_check(m, df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(5.194994994994995, 0.42870445777845795, 'mode')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 30252 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(-10, 10, 1000)\n",
    "pdf = stats.norm.pdf(x, -2, 1) + stats.norm.pdf(x, 5, 0.05) * 0.06\n",
    "mode = pdf.max()\n",
    "mean = pdf.mean()\n",
    "\n",
    "# find the mean of the integrated probablity density function\n",
    "idx_mean = np.argmin(np.abs(np.cumsum(pdf) - np.cumsum(pdf).mean()))\n",
    "idx_mode = pdf.argmax()\n",
    "plt.figure(figsize=(16, 6))\n",
    "plt.plot(x, pdf)\n",
    "plt.vlines([x[idx_mean], x[idx_mode]], 0, mode, lw=0.8, linestyles='--')\n",
    "plt.annotate('mean', (x[idx_mean] + 0.2, mode * 0.9))\n",
    "plt.annotate('mode', (x[idx_mode] + 0.2, mode * 0.9))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "logp = 4,021.6, ||grad|| = 139.89: 100%|██████████| 253/253 [00:00<00:00, 796.77it/s] \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fdd88144198>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 44024 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Find a point estimate of the models parameters\n",
    "with m:\n",
    "    aprox = pm.find_MAP()\n",
    "\n",
    "# Determine g, based on the parameters\n",
    "def det_trend(k, m, delta, t, s, A):\n",
    "    return (k + np.dot(A, delta)) * t + (m + np.dot(A, (-s * delta)))\n",
    "\n",
    "# run function and rescale to original scale\n",
    "g = det_trend(aprox['k'], aprox['m'], aprox['delta'], df['t'], s, A) * df['y'].max()\n",
    "\n",
    "plt.figure(figsize=(16, 6))\n",
    "plt.title('$g(t)$')\n",
    "plt.plot(g)\n",
    "plt.scatter(np.arange(df.shape[0]), df.y, s=0.5, color='black')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def det_trend(k, m, delta, t, s, A):\n",
    "    return (k + np.dot(A, delta)) * t + (m + np.dot(A, (-s * delta)))\n",
    "\n",
    "# k = trace['k'].mean(-1)\n",
    "# m = trace['m'].mean(-1)\n",
    "# delta = trace['delta'].mean(0)\n",
    "\n",
    "k = aprox['k']\n",
    "m = aprox['m']\n",
    "delta = aprox['delta']\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Seasonality\n",
    "\n",
    "$$ s(t) = \\sum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fdd880e2f28>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 45324 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(6)\n",
    "def fourier_series(t, p=365.25, n=10):\n",
    "    # 2 pi n / p\n",
    "    x = 2 * np.pi * np.arange(1, n + 1) / p\n",
    "    # 2 pi n / p * t\n",
    "    x = x * t[:, None]\n",
    "    x = np.concatenate((np.cos(x), np.sin(x)), axis=1)\n",
    "    return x\n",
    "\n",
    "n = 4\n",
    "t = np.arange(1000)\n",
    "beta = np.random.normal(size=2 * n)\n",
    "plt.figure(figsize=(16, 6))\n",
    "plt.plot(fourier_series(t, 365.25, n) @ beta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def seasonality_model(m, df, period='yearly', seasonality_prior_scale=10):\n",
    "    \n",
    "    if period == 'yearly':\n",
    "        n = 10\n",
    "        # rescale the period, as t is also scaled\n",
    "        p = 365.25 / (df['ds'].max() - df['ds'].min()).days\n",
    "    else:  # weekly\n",
    "        n = 3\n",
    "        # rescale the period, as t is also scaled\n",
    "        p = 7 / (df['ds'].max() - df['ds'].min()).days\n",
    "    x = fourier_series(df['t'], p, n)\n",
    "    with m:\n",
    "        beta = pm.Normal(f'beta_{period}', mu=0, sd=seasonality_prior_scale, shape=2 * n)\n",
    "    return x, beta\n",
    "\n",
    "m = pm.Model()\n",
    "\n",
    "with m:\n",
    "    # changepoints_prior_scale is None, so the exponential distribution\n",
    "    # will be used as prior on \\tau.\n",
    "    y, A, s = trend_model(m, df['t'], changepoints_prior_scale=None)\n",
    "    x_yearly, beta_yearly = seasonality_model(m, df, 'yearly')\n",
    "    x_weekly, beta_weekly = seasonality_model(m, df, 'weekly')\n",
    "    \n",
    "    y += det_dot(x_yearly, beta_yearly) + det_dot(x_weekly, beta_weekly)\n",
    "    \n",
    "    sigma = pm.HalfCauchy('sigma', 0.5, testval=1)\n",
    "    obs = pm.Normal('obs', \n",
    "                 mu=y, \n",
    "                 sd=sigma,\n",
    "                 observed=df['y_scaled'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 99724 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sanity_check(m, df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<MultiTrace: 4 chains, 20 iterations, 9 variables>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[<matplotlib.axes._subplots.AxesSubplot object at 0x7fdd78715908>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7fdd78723be0>],\n",
       "       [<matplotlib.axes._subplots.AxesSubplot object at 0x7fdd786b8240>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7fdd786ca860>],\n",
       "       [<matplotlib.axes._subplots.AxesSubplot object at 0x7fdd786d9e80>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7fdd786794e0>],\n",
       "       [<matplotlib.axes._subplots.AxesSubplot object at 0x7fdd7869b898>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7fdd7863ce80>],\n",
       "       [<matplotlib.axes._subplots.AxesSubplot object at 0x7fdd78664518>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7fdd78609b38>],\n",
       "       [<matplotlib.axes._subplots.AxesSubplot object at 0x7fdd785b8198>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7fdd785dd7b8>],\n",
       "       [<matplotlib.axes._subplots.AxesSubplot object at 0x7fdd78581dd8>,\n",
       "        <matplotlib.axes._subplots.AxesSubplot object at 0x7fdd78530438>]],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 605636 characters]\n",
      "text/plain": [
       "<Figure size 864x1008 with 14 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with open('trace.pkl', 'rb') as f:\n",
    "    trace = pickle.load(f)\n",
    "with m:\n",
    "    trace = pm.sample(500)\n",
    "pm.traceplot(trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# with open('trace.pkl', 'wb') as f:\n",
    "#     pickle.dump(trace, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 25)"
      ]
     },
     "execution_count": 139,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.quantile(trace['delta'], [0.025, 0.975], axis=0).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {},
   "outputs": [],
   "source": [
    "# vector distributions\n",
    "beta_yearly = trace['beta_yearly']\n",
    "beta_weekly = trace['beta_weekly']\n",
    "delta = trace['delta']\n",
    "beta_yearly = trace['beta_yearly']\n",
    "# scalar distributions\n",
    "k = trace['k']\n",
    "m = trace['m']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PolyCollection at 0x7fdd78e8c908>"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 223052 characters]\n",
      "text/plain": [
       "<Figure size 1152x1296 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def det_seasonality_posterior(beta, x):\n",
    "    return np.dot(x, beta.T)\n",
    "\n",
    "p = 0.025\n",
    "# vector distributions\n",
    "beta_yearly = trace['beta_yearly']\n",
    "beta_weekly = trace['beta_weekly']\n",
    "delta = trace['delta']\n",
    "\n",
    "# scalar distributions\n",
    "k = trace['k']\n",
    "m = trace['m']\n",
    "\n",
    "# determine the posterior by evaulating all the values in the trace.\n",
    "trend_posterior = ((k + np.dot(A, delta.T)) * df['t'][:, None] + m + np.dot(A, (-s * delta).T)) * df['y'].max()\n",
    "\n",
    "yearly_posterior = det_seasonality_posterior(beta_yearly, x_yearly) * df['y'].max()\n",
    "weekly_posterior = det_seasonality_posterior(beta_weekly, x_weekly) * df['y'].max()\n",
    "\n",
    "date = df['ds'].dt.to_pydatetime()\n",
    "sunday = np.argmax(df['ds'].dt.dayofweek)\n",
    "weekdays = ['sunday', 'monday', 'tuesday', 'wednesday', 'thursday', 'friday', 'saturday']\n",
    "idx_year = np.argmax(df['ds'].dt.dayofyear)\n",
    "\n",
    "plt.figure(figsize=(16, 3*6))\n",
    "b = 411\n",
    "plt.subplot(b)\n",
    "plt.title('total')\n",
    "plt.plot(date, \n",
    "         (trend_posterior + yearly_posterior + weekly_posterior).mean(1), lw=0.5)\n",
    "plt.scatter(date, df['y'], s=0.5, color='black')\n",
    "\n",
    "plt.subplot(b + 1)\n",
    "plt.title('trend')\n",
    "plt.plot(date, trend_posterior.mean(1))\n",
    "quant = np.quantile(trend_posterior, [p, 1 - p], axis=1)\n",
    "plt.fill_between(date, quant[0, :], quant[1, :], alpha=0.25)\n",
    "\n",
    "plt.subplot(b + 2)\n",
    "plt.title('yearly')\n",
    "plt.plot(date[idx_year: idx_year + 365], yearly_posterior.mean(1)[idx_year: idx_year + 365])\n",
    "quant = np.quantile(yearly_posterior, [p, 1 - p], axis=1)\n",
    "plt.fill_between(date[idx_year: idx_year + 365], \n",
    "                 quant[0, idx_year: idx_year + 365], quant[1, idx_year: idx_year + 365], alpha=0.25)\n",
    "\n",
    "plt.subplot(b + 3)\n",
    "plt.title('weekly')\n",
    "plt.plot(weekdays, weekly_posterior.mean(1)[sunday: sunday + 7])\n",
    "quant = np.quantile(weekly_posterior, [p, 1 - p], axis=1)\n",
    "plt.fill_between(weekdays, quant[0, sunday: sunday + 7], \n",
    "                 quant[1, sunday: sunday + 7], alpha=0.25)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.2403221432824845"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace['tau'].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_changepoints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [],
   "source": [
    "days = 356\n",
    "\n",
    "history_points = df.shape[0]\n",
    "probability_changepoint = n_changepoints / history_points\n",
    "\n",
    "future = pd.DataFrame({'ds': pd.date_range(df['ds'].min(), \n",
    "                                           df['ds'].max() + pd.Timedelta(days, 'D'), \n",
    "                        df.shape[0] + days)})\n",
    "future['t'] = (future['ds'] - df['ds'].min()) / (df['ds'].max() - df['ds'].min())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0034423407917383822"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "probability_changepoint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.01706735, 1.0668899 , 1.09815785, 1.10606074])"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "new_changepoints = future['t'][future['t'] > 1].values\n",
    "\n",
    "sample = np.random.random(new_changepoints.shape)\n",
    "new_changepoints = new_changepoints[sample <= probability_changepoint]\n",
    "new_delta = stats.laplace(0, trace['tau'].mean()).rvs(new_changepoints.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100, 350)"
      ]
     },
     "execution_count": 137,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(350, 100)"
      ]
     },
     "execution_count": 151,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_samples = 100\n",
    "new_changepoints = future['t'][future['t'] > 1].values\n",
    "sample = np.random.random((new_changepoints.shape[0], n_samples))\n",
    "(new_changepoints[:, None] * np.ones(n_samples))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.02338951, -1.90641708,  1.77132142,  0.9831247 , -0.31936231,\n",
       "       -1.05786406,  1.55602154, -1.80852354,  0.20379304,  1.1959541 ,\n",
       "       -1.44950353,  0.94848813,  0.5520769 , -2.12981528,  4.31342896,\n",
       "       -4.02723603,  0.11133015,  0.61538896, -0.22614584,  0.57782803,\n",
       "       -0.1589887 ,  0.84737617, -1.25719443,  0.10925113,  0.13380373])"
      ]
     },
     "execution_count": 168,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "delta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(26,)"
      ]
     },
     "execution_count": 181,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(-new_s * new_delta).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 224,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.131345090059634"
      ]
     },
     "execution_count": 224,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.abs(delta).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 228,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2564357806990329"
      ]
     },
     "execution_count": 228,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "trace['tau'].std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 231,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PolyCollection at 0x7fdd6ede9e80>"
      ]
     },
     "execution_count": 231,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 31072 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_samples = 1000\n",
    "days = 150\n",
    "history_points = df.shape[0]\n",
    "probability_changepoint = n_changepoints / history_points\n",
    "\n",
    "future = pd.DataFrame({'ds': pd.date_range(df['ds'].min(), \n",
    "                                           df['ds'].max() + pd.Timedelta(days, 'D'), \n",
    "                        df.shape[0] + days)})\n",
    "future['t'] = (future['ds'] - df['ds'].min()) / (df['ds'].max() - df['ds'].min())\n",
    "\n",
    "# vector distributions\n",
    "beta_yearly = trace['beta_yearly'].mean(0)\n",
    "beta_weekly = trace['beta_weekly'].mean(0)\n",
    "delta = trace['delta'].mean(0)\n",
    "\n",
    "# scalar distributions\n",
    "k = trace['k'].mean()\n",
    "m = trace['m'].mean()\n",
    "\n",
    "trend_forecast = []\n",
    "lambda_ = trace['tau'].mean()\n",
    "\n",
    "for n in range(n_samples):\n",
    "    new_changepoints = future['t'][future['t'] > 1].values\n",
    "    sample = np.random.random(new_changepoints.shape)\n",
    "    new_changepoints = new_changepoints[sample <= probability_changepoint]\n",
    "    new_delta = np.r_[delta, \n",
    "                      stats.laplace(0, lambda_).rvs(new_changepoints.shape[0])]\n",
    "    new_s = np.r_[s, new_changepoints]\n",
    "    new_A = (future['t'][:, None] > new_s) * 1\n",
    "\n",
    "    trend_forecast.append(((k + np.dot(new_A, new_delta)) * future['t']  + (m + np.dot(new_A, (-new_s * new_delta)))) * df['y'].max())\n",
    "    \n",
    "trend_forecast = np.array(trend_forecast)\n",
    "\n",
    "date = future['ds'].dt.to_pydatetime()\n",
    "plt.figure(figsize=(16, 4))\n",
    "plt.title('Trend forecasts uncertainty')\n",
    "plt.plot(date, trend_forecast.mean(0))\n",
    "quant = np.quantile(trend_forecast, [0.025, 0.975], axis=0)\n",
    "plt.fill_between(date, quant[0, :], quant[1, :], alpha=0.25)\n"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
