# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     name: python3
# ---

# + [markdown] id="view-in-github" colab_type="text"
# <a href="https://colab.research.google.com/github/cormach/bayesian_stats_by_b_lambert/blob/master/Chapter_16_q_2.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>

# + id="ws5_orHDrRH6" colab_type="code" colab={}
import os
import json
import shutil
import urllib.request
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# + id="3vnZcMXFrbCB" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 228} outputId="4c7ad99b-b7b4-4b33-f751-bae6726f6b4d"
# Please use the latest version of CmdStanPy
# !pip install --upgrade cmdstanpy

# + id="RXjimxTSr_ER" colab_type="code" colab={}
# Install pre-built CmdStan binary
# (faster than compiling from source via install_cmdstan() function)
tgz_file = 'colab-cmdstan-2.23.0.tar.gz'
tgz_url = 'https://github.com/stan-dev/cmdstan/releases/download/v2.23.0/colab-cmdstan-2.23.0.tar.gz'
if not os.path.exists(tgz_file):
    urllib.request.urlretrieve(tgz_url, tgz_file)
    shutil.unpack_archive(tgz_file)

# + id="lCBvUUUgtaXh" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 35} outputId="d6c0a819-82ed-45a6-ae09-ac336ea417c1"
# Specify CmdStan location via environment variable
os.environ['CMDSTAN'] = './cmdstan-2.23.0'
# Check CmdStan path
from cmdstanpy import CmdStanModel, cmdstan_path
cmdstan_path()

# + id="VdvsYM-t0vtL" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 427} outputId="83d60981-4758-46a9-c08f-53d053a30f19"
# !pip install arviz

# + id="rSHpEIDQ0y3s" colab_type="code" colab={}
import arviz as az

# + [markdown] id="qGmWngiJwWiD" colab_type="text"
# Question 16_2 - hangover

# + id="wPUnG2pXHydF" colab_type="code" colab={}
url='https://raw.githubusercontent.com/alexandrahotti/Solutions-to-A-Students-Guide-to-Bayesian-Statistics-by-Ben-Lambert/master/All_data/stan_hangover.csv'
df = pd.read_csv(url, error_bad_lines=False, header=0)

# + id="fYH_-HnZH7K0" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 295} outputId="067264b4-d5ac-47c8-ca8b-46a85f7e3c3e"
df['date']=pd.to_datetime(df['date'])
fig, ax = plt.subplots()
df.set_index('date').resample('W').sum()['volume'].plot(ax=ax)
df['holiday']=df['holiday']*50
df.set_index('date').resample('W').sum()['holiday'].plot(ax=ax)

# + id="O2HUCkkxTmDy" colab_type="code" colab={}
stan_text = '''data {
    int T;
    real V[T];
    real h[T];
}
parameters {
    real beta0;
    real beta1;
    real<lower=0> sigma;
}
model {
    for(t in 1:T){
      V[t] ~ normal(beta0 + beta1 *h[t], sigma);
    }

    beta0 ~ normal(0, 50);
    beta1 ~ normal(0, 50);
    sigma ~ normal(0, 10);
}
generated quantities {
  real uplift;
  uplift = beta1 / beta0;
}'''
with open('stan_file.stan', 'w') as f:
  f.write(stan_text)

# + id="DbM6Yj_dwIx3" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 391} outputId="a95f358c-45e2-41c6-a995-18b1a23f4e11"
# !cat stan_file.stan

# + id="4_32zKiGmViR" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 66} outputId="bf590d5e-2cf7-4e12-d4fe-3dbc37d3354e"
stan_model = CmdStanModel(stan_file='stan_file.stan')

# + id="Y0t2af4MzGnZ" colab_type="code" colab={}


# + id="yYNxcLa3zOb1" colab_type="code" colab={}
df['holiday']=df['holiday']/50

data = {'V':df.volume.to_numpy(),
        'h':df.holiday.to_numpy(),
        'T':df.shape[0] }

# + id="L91F4UkWzj-i" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 147} outputId="9f2d938a-db62-4c73-c3f4-9d3be7c78383"
stan_posterior=stan_model.sample(data=data)

# + id="-fdhQTBfUncO" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 412} outputId="c9a82d9b-4c1f-4191-aedd-a9474e8bcdf7"
stan_posterior.diagnose()

# + id="BQpVMeYwzniK" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 220} outputId="e0520138-c181-4eae-8204-d912e8f738bc"
stan_posterior.summary().round(decimals=3)

# + id="E-tiKF2ZJP0O" colab_type="code" colab={}
stan_sample = stan_posterior.get_drawset()

# + id="kxYsRoYB1Z_B" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 131} outputId="fd326458-6b78-42a5-d863-ead56cf2ca83"
az_infdata_obj = az.from_cmdstanpy(
    posterior=stan_posterior,
    posterior_predictive="XSim",
    observed_data=data)

az_infdata_obj

# + id="Vy1fyNZR7Zyi" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 895} outputId="68516d46-e5c5-4ca7-c454-de8db5c8921b"
az.plot_autocorr(az_infdata_obj)

# + id="RwUwf-nR8JtZ" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 771} outputId="36d2a801-9b0d-482c-f244-dfce0f391153"
az.plot_pair(az_infdata_obj)

# + id="vvHGMPTS8XrP" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 501} outputId="89b49564-2461-4543-a3ac-a64e7f8638f4"
az.plot_density(az_infdata_obj)

# + id="8mgvuKJs8hjJ" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 747} outputId="d7b3f924-233f-47dc-a131-81be781f8412"
az.plot_trace(az_infdata_obj)

StarcoderdataJupyterScriptsDedupFiltered