贝叶斯深度学习——基于PyMC3的变分推理

%matplotlib inlineimport pymc3 as pmimport theano.tensor as Timport theanoimport sklearnimport numpy as npimport matplotlib.pyplot as pltimport seaborn as sns
sns.set_style('white')from sklearn import datasetsfrom sklearn.preprocessing import scalefrom sklearn.cross_validation import train_test_splitfrom sklearn.datasets import make_moons

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], label='Class 0')
ax.scatter(X[Y==1, 0], X[Y==1, 1], color='r', label='Class 1')
sns.despine(); ax.legend()
ax.set(xlabel='X', ylabel='Y', title='Toy binary classification data set');

# Trick: Turn inputs and outputs into shared variables. 
# It's still the same thing, but we can later change the values of the shared variable # (to switch in the test-data later) and pymc3 will just use the new data. 
# Kind-of like a pointer we can redirect.
# For more info, see: http://deeplearning.net/software/theano/library/compile/shared.html
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)

%%timewith neural_network:
    # Run ADVI which returns posterior means, standard deviations, and the evidence lower bound (ELBO)    v_params = pm.variational.advi(n=50000)Iteration 0 [0%]: ELBO = -368.86Iteration 5000 [10%]: ELBO = -185.65Iteration 10000 [20%]: ELBO = -197.23Iteration 15000 [30%]: ELBO = -203.2Iteration 20000 [40%]: ELBO = -192.46Iteration 25000 [50%]: ELBO = -198.8Iteration 30000 [60%]: ELBO = -183.39Iteration 35000 [70%]: ELBO = -185.04Iteration 40000 [80%]: ELBO = -187.56Iteration 45000 [90%]: ELBO = -192.32Finished [100%]: ELBO = -225.56CPU times: user 36.3 s, sys: 60 ms, total: 36.4 sWall time: 37.2 s

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

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

<matplotlib.text.Text at 0x7fa5dae039b0>

# Replace shared variables with testing setann_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 1pred = ppc['out'].mean(axis=0) > 0.5 

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

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

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

ann_input.set_value(grid_2d)
ann_output.set_value(dummy_out)# Creater posterior predictive samplesppc = pm.sample_ppc(trace, model=neural_network, samples=500)

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])
ax.scatter(X_test[pred==1, 0], X_test[pred==1, 1], color='r')
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])
ax.scatter(X_test[pred==1, 0], X_test[pred==1, 1], color='r')
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)');

# Set back to original data to retrainann_input.set_value(X_train)
ann_output.set_value(Y_train)# Tensors and RV that will be using mini-batchesminibatch_tensors = [ann_input, ann_output]
minibatch_RVs = [out]# Generator that returns mini-batches in each iterationdef 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 = [
    create_minibatch(X_train), 
    create_minibatch(Y_train),
]

total_size = len(Y_train)

%%timewith 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
    )

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

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

pm.traceplot(trace);