notebooks/14_poisreg_with_tfp.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "J33av_vA1Y3m"
   },
   "source": [
    "# Modelling count data with Tensoflow Probability\n",
    "\n",
    "In this notebook you will work with TFP. You will set up regression models that are able to output different conditional probability distributions to model count data. You will define different models with Keras, sklearn and the Tensorflow probability framework and optimize the negative log likelihood (NLL).\n",
    "You compare the performace of the Poisson regression vs. the linear regression on a val dataset.\n",
    "\n",
    "**Dataset:** \n",
    "You work with a camper dataset form https://stats.idre.ucla.edu/r/dae/zip/. The dataset contains data on 250 groups that went to a park. Each group was questioned about how many fish they caught (count), how many children were in the group (child), how many people were in the group (persons), if they used a live bait  and whether or not they brought a camper to the park (camper).\n",
    "You split the data into train and val dataset.\n",
    "\n",
    "**Content:**\n",
    "* Work with different distributions in TFP: Normal and Poisson\n",
    "* Load and split the camper dataset \n",
    "* Fit different regression models to the camper train dataset: linar regression and Poisson regression \n",
    "* Plot the predicted probability distributions (CPD) for two specific datapoints along with their likelihood\n",
    "* Plot the validationdata along with the predicted mean and the 2.5% and 97.5% percentiles of the predicted CPD\n",
    "* Compare the different models based on the val NLL \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "_nIHYY2rSO4I"
   },
   "outputs": [],
   "source": [
    "%tensorflow_version 2.x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "DOhO_5Pt-N9E"
   },
   "source": [
    "#### Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "WLP-37UY1Y31"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from sklearn.model_selection import train_test_split\n",
    "import tensorflow_probability as tfp\n",
    "import tensorflow as tf\n",
    "\n",
    "%matplotlib inline\n",
    "plt.style.use('default')\n",
    "\n",
    "tfd = tfp.distributions\n",
    "tfb = tfp.bijectors\n",
    "print(\"TFP Version\", tfp.__version__)\n",
    "print(\"TF  Version\",tf.__version__)\n",
    "np.random.seed(42)\n",
    "tf.random.set_seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ecQD6xNF1Y38"
   },
   "source": [
    "### Working with a TFP Poisson distribution\n",
    "\n",
    "Here you can see a small example how to work with a Poisson distribution in TFP. The Poisson distribution has only one parameter, often called $\\lambda$ or rate, which defines the mean and the variance of the distribution. We set the rate $\\lambda$ to 2, and plot the probability distribution for the values 0 to 10. Below in the notebook you will define a model to learn this parameter.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "NMzlJ8mO1Y39"
   },
   "outputs": [],
   "source": [
    "dist = tfd.poisson.Poisson(rate = 2) #A\n",
    "vals = np.linspace(0,10,11) #B\n",
    "p = dist.prob(vals) #C\n",
    "print(dist.mean().numpy())  #D\n",
    "print(dist.stddev().numpy())   #E\n",
    "\n",
    "plt.xticks(vals)\n",
    "plt.stem(vals, p,use_line_collection=True)\n",
    "plt.xlabel('Number of Events')\n",
    "plt.ylabel('Probability')\n",
    "plt.show()\n",
    "\n",
    "#A Poisson distribution with parameter rate = 2\n",
    "#B Integer values from 0 to 10 for the x-axis \n",
    "#C Computes the probability for the values\n",
    "#D The mean value yielding 2.0\n",
    "#E The standard deviation yielding sqrt(2.0) = 1.41...\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "18dtcjMG1Y4E"
   },
   "source": [
    "### Loading real count data\n",
    "\n",
    "Here you load the camper data from: https://stats.idre.ucla.edu/r/dae/zip/. The traget variable is the number of fish caught, during a state park visit by a group. You have data of 250 groups that went to the park. Each group was questioned about how many fish they caught (count), how many children were in the group (child), how many people were in the group (persons), if they used a live bait (livebait) and whether or not they brought a camper to the park (camper). This will be the features.\n",
    "You randomly split the data into train and validation dataset (80% train and 20% validation)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "eQuwzZUj1Y4L"
   },
   "outputs": [],
   "source": [
    "# The Fish Data Set\n",
    "# See example 2 from https://stats.idre.ucla.edu/r/dae/zip/ \n",
    "#\"nofish\",\"livebait\",\"camper\",\"persons\",\"child\",\"xb\",\"zg\",\"count\"\n",
    "dat = np.loadtxt('https://raw.githubusercontent.com/tensorchiefs/dl_book/master/data/fish.csv',delimiter=',', skiprows=1)\n",
    "X = dat[...,1:5] #\"livebait\",\"camper\",\"persons\",\"child\"\n",
    "y = dat[...,7]\n",
    "X=np.array(X,dtype=\"float32\")\n",
    "y=np.array(y,dtype=\"float32\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ILaYEVVnBN3F"
   },
   "source": [
    "Let's split the data and look at the counts (how many fish each group caught).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "RCdA9LKL1Y4Y"
   },
   "outputs": [],
   "source": [
    "X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=42,shuffle=True)\n",
    "d = X_train.shape[1]\n",
    "print(X_train.shape, y_train.shape) \n",
    "print(X_val.shape, y_val.shape) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "9mFv-HckTdxq"
   },
   "outputs": [],
   "source": [
    "print(y_val[0:10], y_train[0:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vTcHrKJVByqM"
   },
   "source": [
    "In the following we will look at the number of fish each group caught. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "RVqCbgQZ1Y4e",
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,5))\n",
    "\n",
    "vals, counts = np.unique(y_train, return_counts=True)\n",
    "plt.subplot(1,2,1)\n",
    "plt.stem(vals, counts,use_line_collection=True)\n",
    "plt.xlabel('Count: number of fish caught')\n",
    "plt.ylabel('Frequency')\n",
    "plt.title('Distribution of number of fish caught in training')\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.stem(vals, counts,use_line_collection=True)\n",
    "plt.xlabel('Count: number of fish caught')\n",
    "plt.ylabel('Frequency')\n",
    "plt.xlim(-1,10)\n",
    "plt.title('Zoomed distribution of number of fish caught in training')\n",
    "plt.show()\n",
    "\n",
    "np.max(y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fWdfVAeiPhMP"
   },
   "source": [
    "You see that most of the groups didn't catch any fish at all. Most of the groups were not very successful, but there is one group that was very successful and caught 149 fish!"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "T2LZ1O7SPyZe"
   },
   "source": [
    "Lets pick the two val observations 31 and 33, which you will investigate in the following. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "lNnxCICNwElp"
   },
   "outputs": [],
   "source": [
    "print(X_val[31])#\"livebait\",\"camper\",\"persons\",\"child\"\n",
    "print(X_val[33])#\"livebait\",\"camper\",\"persons\",\"child\"\n",
    "print(y_val[31])#\"number of caught fish\n",
    "print(y_val[33])#\"number of caught fish"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "MAzV_poQQLqe"
   },
   "source": [
    "Group 31 used livebait, had a camper and were 4 persons with one child. They caught 5 fish.  \n",
    "Group 33 used livebait, didn't have a camper and were 4 persons with two childern. They caught 0 fish."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "QeOx8Gk91Y4l"
   },
   "source": [
    "## Linear regression with constant variance\n",
    "\n",
    "In the next few cells you will ignore the fact that you are dealing with count data here and just fit a linear regression model with constant variance to the data. You will fist do this with sklearn and then with keras. You will use the standart MSE loss and calculate the optimal standart deviation to minimize the NLL. Finally, you predict the val data and compare the performance with the RSME, MAE and the NLL. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Y2jVnNH71Y4m"
   },
   "source": [
    "### Linear regression with sklearn \n",
    " \n",
    "Let's fist fit the linear regression with sklean on the training data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "kl72JqqH1Y4o"
   },
   "outputs": [],
   "source": [
    "# The linear regression using non deep learning methods\n",
    "from sklearn.linear_model import LinearRegression\n",
    "model_skl = LinearRegression()\n",
    "res = model_skl.fit(X_train, y_train)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "43AlkxLU1Y4w"
   },
   "source": [
    "In linear regression, we assuming that the $\\sigma$ is constant. To calculate the NLL, we need to estimate this quantity from the training data. The prediction is done on the validation data. Note that we calculate the mean  validation NLL."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "aZrTcLWA1Y4z"
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "# Calculation of the the optimal sigma \n",
    "y_hat_train = model_skl.predict(X_train)\n",
    "n = len(y_hat_train)\n",
    "sigma_hat_2 = np.var(y_train - y_hat_train.flatten(),ddof=2)\n",
    "print('Estimated variance ', sigma_hat_2)\n",
    "print('Estimated standart deviation ', np.sqrt(sigma_hat_2))\n",
    "\n",
    "y_hat = model_skl.predict(X_val) #Prediction on the validationset\n",
    "RMSE_skl = np.sqrt(np.mean((y_val - y_hat.flatten())**2))\n",
    "MAE_skl = np.mean(np.abs(y_val- y_hat.flatten())) \n",
    "\n",
    "NLL_skl =  0.5*np.log(2 * np.pi * sigma_hat_2) + 0.5*np.mean((y_val - y_hat.flatten())**2)/sigma_hat_2\n",
    "print('NLL on training:', 0.5*np.log(2 * np.pi * sigma_hat_2) + 0.5*np.mean((y_train - y_hat_train.flatten())**2)/sigma_hat_2)\n",
    "\n",
    "df1 = pd.DataFrame(\n",
    "          {'RMSE' : RMSE_skl, 'MAE' : MAE_skl, 'NLL (mean)' : NLL_skl}, index=['Linear Regression (sklearn)']\n",
    ")\n",
    "df1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AxwxHcV_1Y4_"
   },
   "source": [
    "### Linear regression with Keras \n",
    " \n",
    "Let's do the same as before with sklearn, but this time you fit a linear regression  model with keras.\n",
    "You have 4 inputs (child , persons livebait, camper) and 1 output (count). Note that you'll use the standart MSE loss."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "J8oUGFcK1Y5C"
   },
   "outputs": [],
   "source": [
    "from tensorflow.keras import Sequential, Model\n",
    "from tensorflow.keras.layers import Dense, Input \n",
    "from tensorflow.keras.optimizers import Adam\n",
    "\n",
    "model_lr = Sequential() \n",
    "model_lr.add(Dense(1,input_dim=d, activation='linear')) \n",
    "model_lr.compile(loss='mean_squared_error',optimizer=Adam(learning_rate=0.01))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "7bZZMADo1Y5K"
   },
   "outputs": [],
   "source": [
    "hist_lr = model_lr.fit(x=X_train, y=y_train, validation_data=(X_val, y_val), epochs=5000, verbose=0, batch_size=len(y_train))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Eo4VNK9J1Y5R"
   },
   "outputs": [],
   "source": [
    "plt.plot(hist_lr.history['loss']) #Note this is the MSE and not the RMSE\n",
    "plt.plot(hist_lr.history['val_loss'])\n",
    "plt.legend(['loss', 'val_loss'])\n",
    "plt.ylabel('MSE')\n",
    "plt.xlabel('Epochs')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "kzvS6_dEV7K2"
   },
   "source": [
    "#### Evaluation of the Performance \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Qyzmr9jLKJiK"
   },
   "outputs": [],
   "source": [
    "# Calculation of the the optimal sigma \n",
    "y_hat_train = model_lr.predict(X_train)\n",
    "sigma_hat_2 = np.var(y_train - y_hat_train.flatten(),ddof=2)\n",
    "print('Estimated variance ', sigma_hat_2)\n",
    "print('Estimated standart deviation ', np.sqrt(sigma_hat_2))\n",
    "\n",
    "y_hat = model_lr.predict(X_val) #Prediction on the validationset\n",
    "RMSE_lr = np.sqrt(np.mean((y_val - y_hat.flatten())**2))\n",
    "MAE_lr = np.mean(np.abs(y_val - y_hat.flatten())) \n",
    "\n",
    "NLL_lr =  0.5*np.log(2 * np.pi * sigma_hat_2) + 0.5*np.mean((y_val - y_hat.flatten())**2)/sigma_hat_2\n",
    "print('NLL on training:', 0.5*np.log(2 * np.pi * sigma_hat_2) + 0.5*np.mean((y_train - y_hat_train.flatten())**2)/sigma_hat_2)\n",
    "\n",
    "df2 = pd.DataFrame(\n",
    "          {'RMSE' : RMSE_lr, 'MAE' : MAE_lr, 'NLL (mean)' : NLL_lr}, index=['Linear Regression (MSE Keras)']\n",
    ")\n",
    "pd.concat([df1,df2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "OTrEASJJ1Y5f"
   },
   "source": [
    "In the pandas dataframe above you see that the RMSE, MAE and the NLL are same. In the next cell you are comparing the coefficients of the keras and sklearn linear regression models. As you can see you get the same results! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "O_KZ6MI81Y5h"
   },
   "outputs": [],
   "source": [
    "print('weights using deep learning:          ',model_lr.get_weights()[0][:,0])\n",
    "print('weights from sklearn:                 ',res.coef_)\n",
    "print('Intercept (bias) using deep learning: ',model_lr.get_weights()[1][0])\n",
    "print('Intercept (bias) using sklearn:       ',res.intercept_)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "scvDQ-mg1Y5r"
   },
   "source": [
    "Let's plot the observed values vs the predicted mean of caught fish on the val dataset. To inicate the CPD you also plot  the 2.5% and 97.5% percentiles of the predicted CPD. You highlight the observations 31 and 33."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "FlzxWOQT1Y5u"
   },
   "outputs": [],
   "source": [
    "y_hat_val=model_lr.predict(X_val)\n",
    "plt.scatter(y_hat_val, y_val,alpha=0.3)\n",
    "plt.scatter(y_hat_val[33], y_val[33],c=\"orange\",marker='o',edgecolors= \"black\")\n",
    "plt.scatter(y_hat_val[31], y_val[31],c=\"orange\",marker='o',edgecolors= \"black\")\n",
    "sort_idx=np.argsort(y_hat_val,axis=0)\n",
    "plt.plot(y_hat_val[sort_idx].flatten(), y_hat_val[sort_idx].flatten()+2*np.sqrt(sigma_hat_2),linestyle='dashed',c=\"black\")\n",
    "plt.plot(y_hat_val[sort_idx].flatten(), y_hat_val[sort_idx].flatten()-2*np.sqrt(sigma_hat_2),linestyle='dashed',c=\"black\")\n",
    "plt.plot(y_hat_val, y_hat_val, c=\"black\")\n",
    "plt.title('Comparison on the validationset')\n",
    "plt.xlabel('predicted average of caught fish')\n",
    "plt.ylabel('observed number of caught fish')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "shnjM81c1pEf"
   },
   "outputs": [],
   "source": [
    "# Let's check the mean of the predicted CPDs for the obeservations nr 31 and 33\n",
    "print(y_hat_val[31])\n",
    "print(y_hat_val[33])\n",
    "# Remember the observed nr of caught fish for the obeservations nr 31 and 33\n",
    "print(y_val[31])\n",
    "print(y_val[33])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Y9R790TuTw-9"
   },
   "source": [
    "Lets check the predicted outcome distribution for the observations 31 and 33."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "bUzQHcEh2fy2"
   },
   "outputs": [],
   "source": [
    "dist = tfd.Normal(loc=y_hat_val,scale=np.sqrt(sigma_hat_2,dtype=\"float32\"))\n",
    "plt.figure(figsize=(14,6))\n",
    "plt.subplot(1,2,1)\n",
    "plt.plot(np.arange(-40,40,0.1),dist.prob(np.arange(-40,40,0.1))[31])\n",
    "plt.vlines(y_hat_val[31], ymin=0, ymax=dist.prob(y_hat_val)[31],linestyle='dashed')\n",
    "plt.vlines(np.expand_dims(y_val,axis=1)[31], ymin=0, ymax=dist.prob(np.expand_dims(y_val,axis=1))[31],linestyle='dotted',color=\"purple\",linewidth=2)\n",
    "plt.xlabel('Number of Events')\n",
    "plt.ylabel('Probability density')\n",
    "plt.title('val observation 31, observed fish=5')\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.plot(np.arange(-40,40,0.1),dist.prob(np.arange(-40,40,0.1))[33])\n",
    "plt.vlines(y_hat_val[33], ymin=0, ymax=dist.prob(y_hat_val)[33],linestyle='dashed')\n",
    "plt.vlines(np.expand_dims(y_val,axis=1)[33], ymin=0, ymax=dist.prob(np.expand_dims(y_val,axis=1))[33],linestyle='dotted',color=\"purple\",linewidth=2)\n",
    "plt.xlabel('Number of Events')\n",
    "plt.ylabel('Probability density')\n",
    "plt.title('val observation 33, observed fish=0')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "UQQlDYPt1Y5_"
   },
   "source": [
    "You can see that the likelihood of the observed values are quite high under the predicted CPDs (dotted line). However, note that the linear model predicts also negative values, which is obviously wrong. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "RC5yQ6Na1Y6u"
   },
   "source": [
    "## Poisson Regression \n",
    "\n",
    "Now you use  the TFP framework and the Poission distribution to model the output of the network as a Poissonian CPD. You will not use any hidden layers in between and the loss will be the NLL. After the fitting, you predict the val data and compare the performance with the linear regression model.\n",
    "$$\n",
    "    Y \\thicksim \\tt{Pois}(exp(w^{T} \\cdot x + b))\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "c1D6mfOl1Y6w"
   },
   "outputs": [],
   "source": [
    "inputs = Input(shape=(X_train.shape[1],))  \n",
    "rate = Dense(1, \n",
    "         activation=tf.exp)(inputs) #A\n",
    "p_y = tfp.layers.DistributionLambda(tfd.Poisson)(rate) #B \n",
    "\n",
    "model_p = Model(inputs=inputs, outputs=p_y) #C\n",
    "\n",
    "\n",
    "def NLL(y_true, y_hat): #D\n",
    "  return -y_hat.log_prob(y_true)\n",
    "\n",
    "model_p.compile(Adam(learning_rate=0.01), loss=NLL)\n",
    "model_p.summary()\n",
    "\n",
    "#A Definition of a single layer with one output\n",
    "#B We use exponential of the output to model the rate\n",
    "#C Glueing the NN and the output layer together. Note that output p_y is a tf.distribution\n",
    "#D The second argument is the output of the model and thus a tfp-distribution. It's as simple as calling log_prob to calculate the log-probability of the observation that is needed to calculate the NLL."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "kheyBbO_1Y6-"
   },
   "outputs": [],
   "source": [
    "hist_p = model_p.fit(x=X_train, y=y_train, validation_data=(X_val, y_val), epochs=2000, verbose=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "kSO3GedG1Y7I",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.plot(hist_p.history['loss'])\n",
    "plt.plot(hist_p.history['val_loss'])\n",
    "plt.legend(['loss', 'val_loss'])\n",
    "plt.xlabel('Epochs')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "3koDxrd5V_KJ"
   },
   "source": [
    "#### Evaluation of the Performance \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "oVjw2n-C1Y7P"
   },
   "outputs": [],
   "source": [
    "#model = Model(inputs=inputs, outputs=p_y.mean()) \n",
    "#y_hat_val = model.predict(X_val).flatten()\n",
    "\n",
    "#New API\n",
    "y_hat_val_dist = model_p(X_val)\n",
    "y_hat_val = y_hat_val_dist.mean().numpy().flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "HeTq4uRalXyN"
   },
   "source": [
    "#### Exercise:\n",
    "Calculate the rmse (root mean square error), the mae (mean absolute error) and the NLL (negative log likelihood) of the Poisson regression model and compare the result to the normal distribution regression model.  \n",
    "What do you observe?  \n",
    "Which model would you choose?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "xHQJXxaTlBev"
   },
   "outputs": [],
   "source": [
    "### your code here \n",
    "\n",
    "rmse=\n",
    "mae= \n",
    "\n",
    "NLL = \n",
    "\n",
    "df3 = pd.DataFrame(\n",
    "         { 'RMSE' : rmse, 'MAE' : mae, 'NLL (mean)' : NLL}, index=['Poisson Regression (TFP)']\n",
    ")\n",
    "pd.concat([df1,df2,df3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vWa4RKhaVVAX"
   },
   "source": [
    "In the pandas dataframe above you see that the RMSE, MAE and the NLL of the diferent models. You see that the Poisson regression outperform the linear regression because of the lower NLL."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "mmdYk1eLWFaN"
   },
   "source": [
    "Let's plot the observed values vs the predicted mean of caught fish on the val dataset. To inicate the CPD you also plot  the 2.5% and 97.5% percentiles of the predicted CPD. You highlight the observations 31 and 33."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ddqrjcSb1Y7T",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from scipy.stats import poisson\n",
    "lower=poisson.ppf(0.025, y_hat_val) #ppf(q, mu, loc=0) Percent point function (inverse of cdf — percentiles) \n",
    "upper=poisson.ppf(0.975, y_hat_val) #ppf(q, mu, loc=0) Percent point function (inverse of cdf — percentiles)\n",
    "\n",
    "#lower=np.percentile(dist.sample(100000),2.5,axis=0)\n",
    "#upper=np.percentile(dist.sample(100000),97.5,axis=0)\n",
    "\n",
    "plt.figure(figsize=(12,6))\n",
    "plt.subplot(1,2,1)\n",
    "plt.scatter(y_hat_val, y_val, alpha=0.3)\n",
    "plt.scatter(y_hat_val[33], y_val[33],c=\"orange\",marker='o',edgecolors= \"black\")\n",
    "plt.scatter(y_hat_val[31], y_val[31],c=\"orange\",marker='o',edgecolors= \"black\")\n",
    "plt.title('Comparison on the validationset')\n",
    "plt.xlabel('predicted average of caught fish')\n",
    "plt.ylabel('observed number of caught fish')\n",
    "plt.plot(y_hat_val[np.argsort(y_hat_val,axis=0)].flatten(), lower[np.argsort(y_hat_val,axis=0)],linestyle='dashed',c=\"black\")\n",
    "plt.plot(y_hat_val[np.argsort(y_hat_val,axis=0)].flatten(), upper[np.argsort(y_hat_val,axis=0)],linestyle='dashed',c=\"black\")\n",
    "plt.plot(y_hat_val, y_hat_val, c=\"black\")\n",
    "\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.scatter(y_hat_val, y_val, alpha=0.3)\n",
    "plt.scatter(y_hat_val[33], y_val[33],c=\"orange\",marker='o',edgecolors= \"black\")\n",
    "plt.scatter(y_hat_val[31], y_val[31],c=\"orange\",marker='o',edgecolors= \"black\")\n",
    "plt.title('Zoomed comparison on the validationset')\n",
    "plt.xlabel('predicted average of caught fish')\n",
    "plt.ylabel('observed number of caught fish')\n",
    "plt.plot(y_hat_val[np.argsort(y_hat_val,axis=0)].flatten(), lower[np.argsort(y_hat_val,axis=0)],linestyle='dashed',c=\"black\")\n",
    "plt.plot(y_hat_val[np.argsort(y_hat_val,axis=0)].flatten(), upper[np.argsort(y_hat_val,axis=0)],linestyle='dashed',c=\"black\")\n",
    "plt.plot(y_hat_val, y_hat_val, c=\"black\")\n",
    "plt.xlim([-0.5,6])\n",
    "plt.ylim([-0.5,6])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "T_5JwyCQWdWm"
   },
   "outputs": [],
   "source": [
    "# Let's check the mean of the predicted CPDs for the obeservations nr 31 and 33\n",
    "print(y_hat_val[31])\n",
    "print(y_hat_val[33])\n",
    "# Remember the observed nr of caught fish for the obeservations nr 31 and 33\n",
    "print(y_val[31])\n",
    "print(y_val[33])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "gHT2pj1kWp6c"
   },
   "source": [
    "Lets check the predicted outcome distribution for the observations 31 and 33."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "3KdtKSlq6CDI"
   },
   "outputs": [],
   "source": [
    "x_ax = np.arange(0,20,1)\n",
    "y_hat_val = y_hat_val_dist.mean().numpy()\n",
    "dist = tfd.Poisson(rate=y_hat_val)\n",
    "probs = dist.prob(value = np.tile(x_ax, (50,1)))\n",
    "\n",
    "plt.figure(figsize=(14,6))\n",
    "plt.subplot(1,2,1)\n",
    "plt.stem(x_ax,probs[31,:],use_line_collection=True)\n",
    "plt.xticks(x_ax)\n",
    "plt.vlines(np.expand_dims(y_val,axis=1)[31], ymin=0, ymax=probs[31,np.int(y_val[31])],linestyle='dotted',color=\"purple\",linewidth=4)\n",
    "\n",
    "plt.xlabel('Number of Events')\n",
    "plt.ylabel('Probability')\n",
    "plt.title('val observation 31, observed fish=5')\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.stem(x_ax,probs[33,:],use_line_collection=True)\n",
    "plt.xticks(x_ax)\n",
    "plt.vlines(np.expand_dims(y_val,axis=1)[33], ymin=0, ymax=probs[33,np.int(y_val[33])],linestyle='dotted',color=\"purple\",linewidth=4)\n",
    "plt.xlabel('Number of Events')\n",
    "plt.ylabel('Probability')\n",
    "plt.title('val observation 33, observed fish=0')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "RP8y4gdz1Y7a"
   },
   "source": [
    "You can see that the likelihood of the observed values are quite high under the predicted CPDs (dotted line). Note that the Poisson CPD does only predict non-negative integer values which is a quite nice property for count data.\n",
    "\n",
    "\n",
    " "
   ]
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "collapsed_sections": [],
   "name": "14_poisreg_with_tfp_new.ipynb",
   "provenance": [],
   "toc_visible": true
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}