I am trying to grasp Bayesain statistics with pymc3
I ran this code for a simple linear regression
#Generating data y=a+bx
import pymc3
import numpy as np
N=1000
alpha,beta, sigma = 2.0, 0.5, 1.0
np.random.seed(47)
X = np.linspace(0, 1, N)
Y = alpha + beta*X + np.random.randn(N)*sigma
#Fitting
linear_model = pymc3.Model()
with linear_model:
alpha = pymc3.Normal('alpha', mu=0, sd=10)
beta = pymc3.Normal('beta', mu=0, sd=10)
sigma = pymc3.HalfNormal('sigma', sd=1)
mu = alpha + beta*X
Y_obs = pymc3.Normal('Y_obs', mu=mu, sd=sigma, observed=Y)
start = pymc3.find_MAP(fmin=optimize.fmin_powell)
step = pymc3.NUTS(scaling=start)
trace = pymc3.sample(500, step, start=start)
I dont understand what does the trace stand for
If I understand the Bayesian theory well enough, there supposed to be a belief
function that gets alpha
,beta
and sigma
and outputs the probability of their combination.
How can I get this belief
structure out of the trace
variables ?
The trace
is the output of the Markov Chain Monte Carlo (MCMC) process. It converges to a distribution (e.g., belief) of your parameters, given the data.
You can view the trace using:
pymc3.traceplot(trace, vars=['alpha', 'beta', 'sigma'])
If you would like to see the individual realizations of your regression along each point of the trace, you can do something like this:
import matplotlib.pyplot as plt
a = trace['alpha']
b = trace['beta']
x = np.linspace(0,1,N)
fig = plt.figure(figsize=(12,4))
ax = fig.add_subplot(1,2,1)
plt.scatter(X,Y, color='g', alpha=0.3)
for i in xrange(500):
y = a[i] + b[i] * x
plt.plot(x, y, 'b', alpha=0.02)
ax = fig.add_subplot(1,2,2)
for i in xrange(500):
y = a[i] + b[i] * x
plt.plot(x, y, 'b', alpha=0.02)
plt.show()
Note: It appears that your code is missing the line:
from scipy import optimize
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With