%matplotlib inline
import numpy as np
import pymc as pm
import matplotlib.pyplot as plt

plt.style.use('fivethirtyeight')

N = 100
x = pm.runiform(-1,1,N)
a_true = 0.5
b_true = 1.2
y_true = np.exp(a_true+b_true*x)
shape_true = 10
y = pm.rgamma(alpha=shape_true,beta=shape_true/y_true,size=N)

plt.scatter(x,y)

# Priors
a = pm.Normal('a', mu=0.0, tau=0.0001)
b = pm.Normal('b', mu=0.0, tau=0.0001)
shape = pm.Uniform('shape', lower=0., upper=100.)

# Model
mu = pm.Lambda('mu', lambda a=a,b=b,shape=shape: np.exp(a+b*x) )

# Likelihood
Yi = pm.Gamma('Yi', alpha=shape, beta=shape/mu, value=y, observed=True)

# Predictive model fit
xnew = np.linspace(min(x),max(x),N)
Ex_mu = pm.Lambda('Ex_mu', lambda a=a,b=b: np.exp(a+b*xnew))

M = pm.MCMC(locals())
M.sample(10000,9000)

plt.hist(M.a.trace())
plt.plot(a_true,1,marker='.',markersize=30)

plt.hist(M.b.trace())
plt.plot(b_true,1,marker='.',markersize=30)

plt.hist(M.shape.trace())
plt.plot(shape_true,1,marker='.',markersize=30)

# Plot data and model
plt.scatter(x,y,c="dodgerblue",s=30)
ymu = np.array([np.median(i) for i in M.Ex_mu.trace().T])
yl95 = np.array([np.percentile(i,2.5) for i in M.Ex_mu.trace().T])
yu95 = np.array([np.percentile(i,97.5) for i in M.Ex_mu.trace().T])
plt.plot(xnew,ymu,c="red",linewidth=1.5)
plt.plot(xnew,yl95,c="grey",linewidth=1)
plt.plot(xnew,yu95,c="grey",linewidth=1)