<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 This example is a Jupyter notebook. You can download it or run it interactively on mybinder.org.
Linear regression
Data
The true parameters of the linear regression
Step1:  Generate data
Step2:  Model
The regressors, that is, the input data
Step3:  Note that we added a column of ones to the regressor matrix for the bias term. We model the slope and the bias term in the same node so we do not factorize between them
Step4:  The first element is the slope which multiplies x and the second element is the bias term which multiplies the constant ones. Now we compute the dot product of X and B
Step5:  The noise parameter
Step6:  The noisy observations
Step7:  Inference
Observe the data
Step8:  Construct the variational Bayesian (VB) inference engine by giving all stochastic nodes
Step9:  Iterate until convergence
Step10:  Results
Create a simple predictive model for new inputs
Step11:  Note that we use the learned node B but create a new regressor array for predictions. Plot the predictive distribution of noiseless function values
Step12:  Note that the above plot shows two standard deviation of the posterior of the noiseless function, thus the data points may lie well outside this range. The red line shows the true linear function. Next, plot the distribution of the noise parameter and the true value, 2−2=0.25
Step13:  The noise level is captured quite well, although the posterior has more mass on larger noise levels (smaller precision parameter values). Finally, plot the distribution of the regression parameters and mark the true value
Step14:  In this case, the true parameters are captured well by the posterior distribution.
Improving accuracy
The model can be improved by not factorizing between B and tau but learning their joint posterior distribution. This requires a slight modification to the model by using GaussianGammaISO node
Step15:  This node contains both the regression parameter vector and the noise parameter. We compute the dot product similarly as before
Step16:  However, Y is constructed as follows
Step17:  Because the noise parameter is already in F_tau we can give a constant one as the second argument. The total noise parameter for Y is the product of the noise parameter in F_tau and one. Now, inference is run similarly as before
<ASSISTANT_TASK:>
Python Code:
import numpy as np
k = 2 # slope
c = 5 # bias
s = 2 # noise standard deviation
# This cell content is hidden from Sphinx-generated documentation
%matplotlib inline
np.random.seed(42)
Explanation: This example is a Jupyter notebook. You can download it or run it interactively on mybinder.org.
Linear regression
Data
The true parameters of the linear regression:
End of explanation
x = np.arange(10)
y = k*x + c + s*np.random.randn(10)
Explanation: Generate data:
End of explanation
X = np.vstack([x, np.ones(len(x))]).T
Explanation: Model
The regressors, that is, the input data:
End of explanation
from bayespy.nodes import GaussianARD
B = GaussianARD(0, 1e-6, shape=(2,))
Explanation: Note that we added a column of ones to the regressor matrix for the bias term. We model the slope and the bias term in the same node so we do not factorize between them:
End of explanation
from bayespy.nodes import SumMultiply
F = SumMultiply('i,i', B, X)
Explanation: The first element is the slope which multiplies x and the second element is the bias term which multiplies the constant ones. Now we compute the dot product of X and B:
End of explanation
from bayespy.nodes import Gamma
tau = Gamma(1e-3, 1e-3)
Explanation: The noise parameter:
End of explanation
Y = GaussianARD(F, tau)
Explanation: The noisy observations:
End of explanation
Y.observe(y)
Explanation: Inference
Observe the data:
End of explanation
from bayespy.inference import VB
Q = VB(Y, B, tau)
Explanation: Construct the variational Bayesian (VB) inference engine by giving all stochastic nodes:
End of explanation
Q.update(repeat=1000)
Explanation: Iterate until convergence:
End of explanation
xh = np.linspace(-5, 15, 100)
Xh = np.vstack([xh, np.ones(len(xh))]).T
Fh = SumMultiply('i,i', B, Xh)
Explanation: Results
Create a simple predictive model for new inputs:
End of explanation
import bayespy.plot as bpplt
bpplt.pyplot.figure()
bpplt.plot(Fh, x=xh, scale=2)
bpplt.plot(y, x=x, color='r', marker='x', linestyle='None')
bpplt.plot(k*xh+c, x=xh, color='r');
Explanation: Note that we use the learned node B but create a new regressor array for predictions. Plot the predictive distribution of noiseless function values:
End of explanation
bpplt.pyplot.figure()
bpplt.pdf(tau, np.linspace(1e-6,1,100), color='k')
bpplt.pyplot.axvline(s**(-2), color='r');
Explanation: Note that the above plot shows two standard deviation of the posterior of the noiseless function, thus the data points may lie well outside this range. The red line shows the true linear function. Next, plot the distribution of the noise parameter and the true value, 2−2=0.25:
End of explanation
bpplt.pyplot.figure();
bpplt.contour(B, np.linspace(1,3,1000), np.linspace(1,9,1000),
              n=10, colors='k');
bpplt.plot(c, x=k, color='r', marker='x', linestyle='None',
           markersize=10, markeredgewidth=2)
bpplt.pyplot.xlabel(r'$k$');
bpplt.pyplot.ylabel(r'$c$');
Explanation: The noise level is captured quite well, although the posterior has more mass on larger noise levels (smaller precision parameter values). Finally, plot the distribution of the regression parameters and mark the true value:
End of explanation
from bayespy.nodes import GaussianGamma
B_tau = GaussianGamma(np.zeros(2), 1e-6*np.identity(2), 1e-3, 1e-3)
Explanation: In this case, the true parameters are captured well by the posterior distribution.
Improving accuracy
The model can be improved by not factorizing between B and tau but learning their joint posterior distribution. This requires a slight modification to the model by using GaussianGammaISO node:
End of explanation
F_tau = SumMultiply('i,i', B_tau, X)
Explanation: This node contains both the regression parameter vector and the noise parameter. We compute the dot product similarly as before:
End of explanation
Y = GaussianARD(F_tau, 1)
Explanation: However, Y is constructed as follows:
End of explanation
Y.observe(y)
Q = VB(Y, B_tau)
Q.update(repeat=1000)
Explanation: Because the noise parameter is already in F_tau we can give a constant one as the second argument. The total noise parameter for Y is the product of the noise parameter in F_tau and one. Now, inference is run similarly as before:
End of explanation
<END_TASK>