In [1]:
import pystan
import numpy as np
import pandas as pd
import os
import pickle
import matplotlib.pyplot as plt
import tensorflow as tf
import edward as ed
In [2]:
%matplotlib inline
In [3]:
model_code = """
data {
int N; //the number of observations
int N2; //the size of the new_X matrix
int K; //the number of columns in the model matrix
vector[N] y; //the response
matrix[N,K] X; //the model matrix
matrix[N2,K] new_X; //the matrix for the predicted values
}
parameters {
vector[K] beta; //the regression parameters
real sigma; //the standard deviation
}
transformed parameters {
vector[N] linpred;
linpred = X*beta;
}
model {
beta[1] ~ cauchy(0,10); //prior for the intercept following Gelman 2008
for(i in 2:K)
beta[i] ~ cauchy(0,2.5);//prior for the slopes following Gelman 2008
y ~ normal(linpred,sigma);
}
generated quantities {
vector[N2] y_pred;
y_pred = new_X*beta; //the y values predicted by the model
}
"""
In [4]:
model_name = 'linear_regression'
pkl_file = model_name + '.pkl'
if os.path.isfile(pkl_file):
# Reuse previously compiled model
sm = pickle.load(open(pkl_file, 'rb'))
else:
# Compile and sample model
sm = pystan.StanModel(model_code=model_code, model_name=model_name)
with open(pkl_file, 'wb') as f:
pickle.dump(sm, f)
In [5]:
# Test data
N = 20
beta = np.array([0.2, 1.5]).reshape((-1,1))
x = np.linspace(0, 1, N).reshape((-1, 1))
X = np.column_stack([np.ones(x.shape), x])
y = X.dot(beta) + np.random.normal(0, 0.1, size=N).reshape((-1, 1))
In [6]:
plt.plot(x, y, marker='.', linestyle='none')
Out[6]:
In [7]:
# Plot fit
N2 = 2
new_x = np.linspace(0, 1, N2).reshape((-1, 1))
new_X = np.column_stack([np.ones(N2).reshape((-1, 1)), new_x])
In [8]:
data = {
'N': N,
'X': X,
'K': X.shape[1],
'y': y.flatten(),
'N2': N2,
'new_X': new_X
}
In [9]:
res = sm.optimizing(data=data)
In [10]:
# Plot MLE fit
plt.plot(x, y, marker='.', linestyle='none')
plt.plot(new_x, res['y_pred'])
plt.show()
In [11]:
res_vb = sm.vb(data=data, pars=['beta', 'y_pred'])
out_file = res_vb['args']['sample_file']
df = pd.read_csv(out_file, comment='#').iloc[2:, :]
In [12]:
df.describe()
Out[12]:
In [13]:
df.head()
Out[13]:
In [14]:
df[['beta.1', 'beta.2']].hist(bins=50)
Out[14]:
In [15]:
# Plot with generated values after Variational Bayes fit
plt.plot(x, y, marker='.', linestyle='none')
plt.plot(np.tile(new_x, (1, len(df))), df[['y_pred.1', 'y_pred.2']].T.values, alpha=0.005, c='k')
plt.show()
Try to replicate Bayesian linear regression using Edward - see http://edwardlib.org/tutorials/supervised-regression
In [16]:
Beta_ed = ed.models.Normal(loc=tf.zeros([2, 1]), scale=tf.ones([2, 1]))
X_ed = tf.placeholder(tf.float32, [N, 2]) # X.astype(np.float32)
y_ed = ed.models.Normal(loc=tf.matmul(X_ed, Beta_ed), scale=0.5)
In [17]:
# Different prior - doesn't seem to be a Cauchy distribution in Edward
qBeta = ed.models.Normal(
loc=tf.Variable(tf.zeros([2, 1])),
scale=tf.nn.softplus(tf.Variable(tf.zeros([2, 1]))))
In [18]:
inference = ed.KLqp({Beta_ed: qBeta}, data={X_ed: X, y_ed: y})
In [19]:
inference.run()
In [20]:
n_sample = 1000
df_ed = pd.DataFrame(qBeta.sample(n_sample).eval().reshape((n_sample, -1)), columns=['beta.1', 'beta.2'])
In [21]:
df_ed[['beta.1', 'beta.2']].hist(bins=50)
plt.show()
In [22]:
# Plot with generated values after Variational Bayes fit
plt.plot(x, y, marker='.', linestyle='none')
plt.plot(np.tile(new_x, (1, len(df_ed))),
new_X.dot(df_ed[['beta.1', 'beta.2']].values.T), alpha=0.005, c='k')
plt.show()
In [ ]: