PyMC3 - 贝叶斯神经网络

当前机器学习趋势

大规模概率编程

深度学习

架起深度学习和概率编程的桥梁

PyMC3中的贝叶斯深网络

生成数据

%matplotlib inline
import pymc3 as pm
import theano.tensor as T
import theano
import sklearn
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('white')
from sklearn import datasets
from sklearn.preprocessing import scale
from sklearn.cross_validation import train_test_split
from sklearn.datasets import make_moons

/Users/xujie/.pyenv/versions/anaconda3-4.3.1/lib/python3.6/site-packages/sklearn/cross_validation.py:44: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. Also note that the interface of the new CV iterators are different from that of this module. This module will be removed in 0.20.
"This module will be removed in 0.20.", DeprecationWarning)

X, Y = make_moons(noise=0.2, random_state=0, n_samples=1000)
X = scale(X)
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=.5)

fig, ax = plt.subplots()
ax.scatter(X[Y==0, 0], X[Y==0, 1], color='b',edgecolors='k',alpha=0.6, label='Class 0')
ax.scatter(X[Y==1, 0], X[Y==1, 1], color='r',edgecolors='k',alpha=0.6, label='Class 1')
sns.despine(); ax.legend()
ax.set(xlabel='X', ylabel='Y', title='Toy binary classification data set');

模型定义

ann_input = theano.shared(X_train)
ann_output = theano.shared(Y_train)

n_hidden = 5

# Initialize random weights between each layer
init_1 = np.random.randn(X.shape[1], n_hidden)
init_2 = np.random.randn(n_hidden, n_hidden)
init_out = np.random.randn(n_hidden)

with pm.Model() as neural_network:
# Weights from input to hidden layer
weights_in_1 = pm.Normal('w_in_1', 0, sd=1,
shape=(X.shape[1], n_hidden),
testval=init_1)

# Weights from 1st to 2nd layer
weights_1_2 = pm.Normal('w_1_2', 0, sd=1,
shape=(n_hidden, n_hidden),
testval=init_2)

# Weights from hidden layer to output
weights_2_out = pm.Normal('w_2_out', 0, sd=1,
shape=(n_hidden,),
testval=init_out)

# Build neural-network using tanh activation function
act_1 = T.tanh(T.dot(ann_input, weights_in_1))
act_2 = T.tanh(T.dot(act_1, weights_1_2))
act_out = T.nnet.sigmoid(T.dot(act_2, weights_2_out))

# Binary classification -> Bernoulli likelihood
out = pm.Bernoulli('out',
act_out,
observed=ann_output)

/Users/xujie/.pyenv/versions/anaconda3-4.3.1/lib/python3.6/site-packages/theano/tensor/basic.py:2146: UserWarning: theano.tensor.round() changed its default from `half_away_from_zero` to `half_to_even` to have the same default as NumPy. Use the Theano flag `warn.round=False` to disable this warning.
"theano.tensor.round() changed its default from"

变分推理：缩放模型的复杂性

NUTS

%%time

with neural_network:
# Run ADVI which returns posterior means, standard deviations, and the evidence lower bound (ELBO)
v_params = pm.variational.advi(n=50000)

Average ELBO = -151.63: 100%|██████████| 50000/50000 [00:13<00:00, 3790.27it/s]
Finished [100%]: Average ELBO = -136.89

CPU times: user 15.5 s, sys: 695 ms, total: 16.1 s
Wall time: 18 s

plt.plot(v_params.elbo_vals);
plt.ylabel('ELBO');
plt.xlabel('iteration');

sample_vp()

with neural_network:
trace = pm.variational.sample_vp(v_params, draws=5000)

100%|██████████| 5000/5000 [00:00<00:00, 10687.44it/s]

sample_ppc()

# Replace shared variables with testing set
ann_input.set_value(X_test)
ann_output.set_value(Y_test)

# Creater posterior predictive samples
ppc = pm.sample_ppc(trace, model=neural_network, samples=500)

# Use probability of > 0.5 to assume prediction of class 1
pred = ppc['out'].mean(axis=0) > 0.5

100%|██████████| 500/500 [00:03<00:00, 139.67it/s]

fig, ax = plt.subplots()
ax.scatter(X_test[pred==0, 0], X_test[pred==0, 1], color='b',edgecolors='k',alpha=0.6)
ax.scatter(X_test[pred==1, 0], X_test[pred==1, 1], color='r',edgecolors='k',alpha=0.6)
sns.despine()
ax.set(title='Predicted labels in testing set', xlabel='X', ylabel='Y');

print('Accuracy = {}%'.format((Y_test == pred).mean() * 100))

Accuracy = 95.6%

分类器分析

grid = np.mgrid[-3:3:100j,-3:3:100j]
grid_2d = grid.reshape(2, -1).T
dummy_out = np.ones(grid.shape[1], dtype=np.int8)

ann_input.set_value(grid_2d)
ann_output.set_value(dummy_out)

# 对后验估计进行采样
ppc = pm.sample_ppc(trace, model=neural_network, samples=500)

100%|██████████| 500/500 [00:04<00:00, 101.43it/s]

绘制预测平面如下

cmap = sns.diverging_palette(250, 12, s=85, l=25, as_cmap=True)
fig, ax = plt.subplots(figsize=(10, 6))
contour = ax.contourf(*grid, ppc['out'].mean(axis=0).reshape(100, 100), cmap=cmap)
ax.scatter(X_test[pred==0, 0], X_test[pred==0, 1], color='b', edgecolors='k',alpha=0.6)
ax.scatter(X_test[pred==1, 0], X_test[pred==1, 1], color='r', edgecolors='k',alpha=0.6)
cbar = plt.colorbar(contour, ax=ax)
_ = ax.set(xlim=(-3, 3), ylim=(-3, 3), xlabel='X', ylabel='Y');
cbar.ax.set_ylabel('Posterior predictive mean probability of class label = 0');

预测中的不确定性

cmap = sns.cubehelix_palette(light=1, as_cmap=True)
fig, ax = plt.subplots(figsize=(10, 6))
contour = ax.contourf(*grid, ppc['out'].std(axis=0).reshape(100, 100), cmap=cmap)
ax.scatter(X_test[pred==0, 0], X_test[pred==0, 1], color='b', edgecolors='k',alpha=0.6)
ax.scatter(X_test[pred==1, 0], X_test[pred==1, 1], color='r', edgecolors='k',alpha=0.6)
cbar = plt.colorbar(contour, ax=ax)
_ = ax.set(xlim=(-3, 3), ylim=(-3, 3), xlabel='X', ylabel='Y');
cbar.ax.set_ylabel('Uncertainty (posterior predictive standard deviation)');

微批AVDI（Mini-batch AVDI）：可扩展数据大小

from six.moves import zip

# Set back to original data to retrain
ann_input.set_value(X_train)
ann_output.set_value(Y_train)

# Tensors and RV that will be using mini-batches
minibatch_tensors = [ann_input, ann_output]
minibatch_RVs = [out]

# Generator that returns mini-batches in each iteration
def create_minibatch(data):
rng = np.random.RandomState(0)

while True:
# Return random data samples of set size 100 each iteration
ixs = rng.randint(len(data), size=50)
yield data[ixs]

minibatches = zip(
create_minibatch(X_train),
create_minibatch(Y_train),
)

total_size = len(Y_train)

%%time

with neural_network:
# Run advi_minibatch
v_params = pm.variational.advi_minibatch(
n=50000, minibatch_tensors=minibatch_tensors,
minibatch_RVs=minibatch_RVs, minibatches=minibatches,
total_size=total_size, learning_rate=1e-2, epsilon=1.0
)

Average ELBO = -123.12: 100%|██████████| 50000/50000 [00:16<00:00, 3043.69it/s]
Finished minibatch ADVI: ELBO = -121.76

CPU times: user 20.5 s, sys: 371 ms, total: 20.9 s
Wall time: 23.9 s

with neural_network:
trace = pm.variational.sample_vp(v_params, draws=5000)

100%|██████████| 5000/5000 [00:00<00:00, 9914.64it/s]

pm.traceplot(trace);