<jupyter_start><jupyter_text>*This example is a Jupyter notebook. You can download it or run it interactively on mybinder.org.* Linear regression DataThe true parameters of the linear regression:<jupyter_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)<jupyter_output><empty_output><jupyter_text>Generate data:<jupyter_code>x = np.arange(10)
y = k*x + c + s*np.random.randn(10)<jupyter_output><empty_output><jupyter_text>ModelThe regressors, that is, the input data:<jupyter_code>X = np.vstack([x, np.ones(len(x))]).T<jupyter_output><empty_output><jupyter_text>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:<jupyter_code>from bayespy.nodes import GaussianARD
B = GaussianARD(0, 1e-6, shape=(2,))<jupyter_output><empty_output><jupyter_text>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:<jupyter_code>from bayespy.nodes import SumMultiply
F = SumMultiply('i,i', B, X)<jupyter_output><empty_output><jupyter_text>The noise parameter:<jupyter_code>from bayespy.nodes import Gamma
tau = Gamma(1e-3, 1e-3)<jupyter_output><empty_output><jupyter_text>The noisy observations:<jupyter_code>Y = GaussianARD(F, tau)<jupyter_output><empty_output><jupyter_text>InferenceObserve the data:<jupyter_code>Y.observe(y)<jupyter_output><empty_output><jupyter_text>Construct the variational Bayesian (VB) inference engine by giving all stochastic nodes:<jupyter_code>from bayespy.inference import VB
Q = VB(Y, B, tau)<jupyter_output><empty_output><jupyter_text>Iterate until convergence:<jupyter_code>Q.update(repeat=1000)<jupyter_output>Iteration 1: loglike=-4.080970e+01 (0.003 seconds)
Iteration 2: loglike=-4.057132e+01 (0.003 seconds)
Iteration 3: loglike=-4.056572e+01 (0.003 seconds)
Iteration 4: loglike=-4.056551e+01 (0.003 seconds)
Converged at iteration 4.<jupyter_text>ResultsCreate a simple predictive model for new inputs:<jupyter_code>xh = np.linspace(-5, 15, 100)
Xh = np.vstack([xh, np.ones(len(xh))]).T
Fh = SumMultiply('i,i', B, Xh)<jupyter_output><empty_output><jupyter_text>Note that we use the learned node B but create a new regressor array for predictions. Plot the predictive distribution of noiseless function values:<jupyter_code>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');<jupyter_output><empty_output><jupyter_text>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:<jupyter_code>bpplt.pyplot.figure()
bpplt.pdf(tau, np.linspace(1e-6,1,100), color='k')
bpplt.pyplot.axvline(s**(-2), color='r');<jupyter_output><empty_output><jupyter_text>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:<jupyter_code>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$');<jupyter_output><empty_output><jupyter_text>In this case, the true parameters are captured well by the posterior distribution. Improving accuracyThe 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:<jupyter_code>from bayespy.nodes import GaussianGamma
B_tau = GaussianGamma(np.zeros(2), 1e-6*np.identity(2), 1e-3, 1e-3)<jupyter_output><empty_output><jupyter_text>This node contains both the regression parameter vector and the noise parameter. We compute the dot product similarly as before:<jupyter_code>F_tau = SumMultiply('i,i', B_tau, X)<jupyter_output><empty_output><jupyter_text>However, Y is constructed as follows:<jupyter_code>Y = GaussianARD(F_tau, 1)<jupyter_output><empty_output><jupyter_text>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:<jupyter_code>Y.observe(y)
Q = VB(Y, B_tau)
Q.update(repeat=1000)<jupyter_output>Iteration 1: loglike=-4.130423e+01 (0.002 seconds)
Iteration 2: loglike=-4.130423e+01 (0.002 seconds)
Converged at iteration 2.
{'fraction_non_alphanumeric': 0.10689931809065382, 'fraction_numerical': 0.031287605294825514, 'mean_word_length': 6.915873015873016, 'pattern_counts': {'":': 0, '<': 71, '<?xml version=': 0, '>': 71, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '57194', 'license': 'MIT', 'n_tokens_mistral': 1824, 'n_tokens_neox': 1697, 'n_words': 491, 'path': 'doc/source/examples/regression.ipynb', 'repo_name': 'jeevan449/NewDoc', 'size': '4986'}