deep-regression.ipynb

Kernel: bayesian

import pymc3 as pm
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import theano.tensor as tt
import theano

%load_ext autoreload
%autoreload 2
%matplotlib inline
%config InlineBackend.figure_format = 'retina'

df = pd.read_csv("../datasets/bikes/hour.csv")
df

feature_cols = ["workingday", "holiday", "temp", "atemp", "hum", "windspeed"]
out_col = ["cnt"]
df[out_col]

X = pm.floatX(df[feature_cols])
Y = pm.floatX(df[out_col].apply(np.log10))
n_hidden = X.shape[1]

with pm.Model() as nn_model:
    w1 = pm.Normal("w1", mu=0, sd=1, shape=(X.shape[1], n_hidden))
    w2 = pm.Normal("w2", mu=0, sd=1, shape=(n_hidden, 1))

    b1 = pm.Normal("b1", mu=0, sd=1, shape=(n_hidden,))
    b2 = pm.Normal("b2", mu=0, sd=1, shape=(1,))

    a1 = pm.Deterministic("a1", tt.nnet.relu(tt.dot(X, w1) + b1))
    a2 = pm.Deterministic("a2", tt.dot(a1, w2) + b2)

    output = pm.Normal("likelihood", mu=a2, observed=Y)

with pm.Model() as three_layer_model:
    w1 = pm.Normal("w1", mu=0, sd=1, shape=(X.shape[1], n_hidden))
    w2 = pm.Normal("w2", mu=0, sd=1, shape=(n_hidden, n_hidden))
    w3 = pm.Normal("w3", mu=0, sd=1, shape=(n_hidden, 1))

    b1 = pm.Normal("b1", mu=0, sd=1, shape=(n_hidden,))
    b2 = pm.Normal("b2", mu=0, sd=1, shape=(n_hidden,))
    b3 = pm.Normal("b3", mu=0, sd=1, shape=(1,))

    a1 = pm.Deterministic("a1", tt.nnet.relu(tt.dot(X, w1) + b1))
    a2 = pm.Deterministic("a2", tt.nnet.relu(tt.dot(a1, w2) + b2))
    a3 = pm.Deterministic("a3", tt.dot(a2, w3) + b3)

    sd = pm.HalfCauchy("sd", beta=1)

    output = pm.Normal("likelihood", mu=a3, sd=sd, observed=Y)

with pm.Model() as linreg_model:
    w1 = pm.Normal("w1", mu=0, sd=1, shape=(X.shape[1], 1))
    b1 = pm.Normal("b1", mu=0, sd=1, shape=(1,))
    a1 = pm.Deterministic("a1", tt.dot(X, w1) + b1)

    sd = pm.HalfCauchy("sd", beta=1)

    output = pm.Normal("likelihood", mu=a1, sd=sd, observed=Y)

with linreg_model:
    s = theano.shared(pm.floatX(1.1))
    inference = pm.ADVI(cost_part_grad_scale=s, learning_rate=0.01)
    approx = pm.fit(200000, method=inference)

plt.plot(inference.hist)

with linreg_model:
    trace = approx.sample(2000)

pm.traceplot(trace, varnames=["w1", "b1"])

with linreg_model:
    samps = pm.sample_ppc(trace)

samps["likelihood"].std(axis=0)

samps["likelihood"].mean(axis=0)

from sklearn.metrics import mean_squared_error as mse

mse(Y, samps["likelihood"].mean(axis=0))

plt.scatter(samps["likelihood"].mean(axis=0).squeeze(), Y.values)