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

# +
import numpy as np
import pandas as pd
import pymc3 as pm
import theano.tensor as tt
import theano
from scipy.stats import norm

import seaborn as sns
import matplotlib.pyplot as plt
# -

# ### Generate Data

# +
n_days = 1000
gpd = 1

true_logλ = 1.5
true_Δlogλ_μ = 0.001
true_Δlogλ_σ = 0.005

lv_rows = []
data_rows = []

for d in range(0, n_days):
    true_Δlogλ = np.random.normal(true_Δlogλ_μ, true_Δlogλ_σ)
    true_logλ = true_logλ + true_Δlogλ
    lv_rows.append([d, true_logλ, true_Δlogλ])
    for _ in range(0,gpd):
        s = np.random.poisson(np.exp(true_logλ))
        data_rows.append([d, s])
# -

lv = pd.DataFrame(lv_rows, columns=['day', 'logλ', 'Δlogλ'])
data = pd.DataFrame(data_rows, columns=['day', 's'])

lv.head()

data.head()

plt.plot(lv['day'], lv['logλ'], color='skyblue', linewidth=1, label='true logλ')
plt.legend()


def get_posteriors(trace):
    p = {}
    logλ_μ, logλ_σ = norm.fit(trace['logλ'])
    p['logλ'] = [logλ_μ, logλ_σ]
    Δlogλ_μ, Δlogλ_σ = norm.fit(trace['Δlogλ'])
    p['Δlogλ'] = [0.0, Δlogλ_σ]
    return p


# +
iv_rows = []

p = {
    'logλ': [1.45, 0.1],
    'Δlogλ': [0.0, 0.005]
}
    
for d in range(0, 500):
    s_obs = data[data['day'] == d]['s']

    with pm.Model() as model:
        logλ_init = pm.Normal('logλ_init', mu=p['logλ'][0], sd=p['logλ'][1])
        Δlogλ = pm.Normal('Δlogλ', mu=p['Δlogλ'][0], sd=p['Δlogλ'][1])
        logλ = pm.Deterministic('logλ', logλ_init + Δlogλ)

        # Likelihood of observed data
        s = pm.Poisson('s', mu=np.exp(logλ), observed=s_obs)

        trace = pm.sample(1000, tune=1000, cores=3)

        p = get_posteriors(trace)
        logλ_μ, logλ_σ = p['logλ']
        Δlogλ_μ, Δlogλ_σ = p['Δlogλ']
        iv_rows.append([d, logλ_μ, logλ_σ, Δlogλ_μ, Δlogλ_σ])
# -

iv = pd.DataFrame(iv_rows, columns=['day', 'logλ_μ', 'logλ_σ', 'Δlogλ_μ', 'Δlogλ_σ'])

iv.tail()

plt.plot(lv.iloc[0:500,:]['day'], lv.iloc[0:500,:]['logλ'], color='skyblue', linewidth=1, label='true logλ')
plt.plot(iv['day'], iv['logλ_μ'], color='red', linewidth=1, label='inferred logλ')
plt.legend()



StarcoderdataJupyterScriptsDedupFiltered