%matplotlib inline
from IPython.display import Image
from IPython.html.widgets import interact

Image('Lighthouse_schematic.jpg',width=500)

xi = linspace(-10,10,150)
alpha = 1
f = lambda x: 1/(pi*(1+(x-alpha)**2))
plot(xi,f(xi))
xlabel('$x$',fontsize=24)
ylabel(r'$f_{\alpha}(x)$',fontsize=24);
vlines(alpha,0,.35,linestyles='--',lw=4.)
grid()

beta  =alpha = 1
theta_samples=((2*np.random.rand(250)-1)*pi/2)
x_samples = alpha+beta*np.tan(theta_samples)
hist(x_samples);

import sympy as S
a=S.symbols('alpha',real=True)
x = S.symbols('x',real=True)
# form derivative of the log-likelihood
dL=sum((xk-a)/(1+(xk-a)**2)  for xk in x_samples)
S.plot(dL,(a,-15,15),xlabel=r'$\alpha$',ylabel=r'$dL(\alpha)/d\alpha$')

from scipy import optimize
from scipy.optimize import fmin_l_bfgs_b

# convert sympy function to numpy version with lambdify
alpha_x  = fmin_l_bfgs_b(S.lambdify(a,(dL)**2),0,bounds=[(-3,3)],approx_grad=True)

print alpha_x

print 'alpha using average =',x_samples.mean()
print 'maximum likelihood estimate = ', alpha_x[0]

def run_trials(n=100):
    o=[]
    for i in range(100):
        theta_samples=((2*np.random.rand(250)-1)*pi/2)
        x_samples = alpha+beta*np.tan(theta_samples)
        o.append(x_samples)
    return np.array(o)

o= run_trials()

hist(o[np.where(abs(o)<200)]);

plot(range(100,10000,100),[o[np.where(abs(o)<i)].mean() for i in range(100,10000,100)],'-o')
xlabel('width of neighborhood around zero')
ylabel('value of average estimate')

def ML_estimator(x_samples):
    dL=sum((xk-a)/(1+(xk-a)**2)  for xk in x_samples)
    a_x  = fmin_l_bfgs_b(S.lambdify(a,(dL)**2),0,bounds=[(-3,3)],approx_grad=True)
    return a_x[0]

# run maximum likelihood estimator on synthetic data we generated earlier
# Beware this may take a long time!
v= np.hstack([ML_estimator(o[i,:]) for i in range(o.shape[0])])

vmed= np.hstack([np.median(o[i,:]) for i in range(o.shape[0])])
vavg =  np.hstack([np.mean(o[i,:]) for i in range(o.shape[0])])
fig,ax = subplots()
ax.plot(v,'-o',label='ML')
ax.plot(vmed,'gs-',label='median')
ax.plot(vavg,'r^-',label='mean')
ax.axis(ymax=2,ymin=0)
ax.legend(loc=(1,0))
ax.set_xlabel('Trial Index')
ax.set_ylabel('Value of Estimator')

fig,ax = subplots()
ii= np.argsort(v)
ax.plot(v[ii],vmed[ii],'o',alpha=.3)
axs = ax.axis()
ax.plot(linspace(0,2,10),linspace(0,2,10))
ax.axis(axs)
ax.set_aspect(1)
ax.set_xlabel('ML estimate')
ax.set_ylabel('Median estimate')

fig,ax = subplots()
ax.hist(v,10,alpha=.3,label='ML')
ax.hist(vmed,10,alpha=.3,label='median')
ax.legend(loc=(1,0))

alpha  = a
alphamx,alphamn=3,-3
g = f(x_samples[0])
xi = linspace(alphamn,alphamx,100)
mxval = S.lambdify(a,g)(xi).max()
plot(xi,S.lambdify(a,g)(xi),x_samples[0],mxval*1.1,'o')

xi = linspace(alphamn,alphamx,100)
palpha=S.Piecewise((1,abs(x)<3),(0,True))

def slide_figure(n=0):
    fig,ax=subplots()
    palpha=S.Piecewise((1,abs(x)<3),(0,True))
    if n==0:
        ax.plot(xi,[S.lambdify(x,palpha)(i) for i in xi])
        ax.plot(x_samples[0],1.1,'o',ms=10.)
    else:
        g = S.prod(map(f,x_samples[:n]))
        mxval = S.lambdify(a,g)(xi).max()*1.1
        ax.plot(xi,S.lambdify(a,g)(xi))
        ax.plot(x_samples[:n],mxval*ones(n),'o',color='gray',alpha=0.3)
        ax.plot(x_samples[n],mxval,'o',ms=10.)
        ax.axis(ymax=mxval*1.1)
    ax.set_xlabel(r'$\alpha$',fontsize=18)
    ax.axis(xmin=-17,xmax=17)

interact(slide_figure, n=(0,15,1));

fig,axs = subplots(2,6,sharex=True)
fig.set_size_inches((10,3))
for n,ax in enumerate(axs.flatten()):
    if n==0:
        ax.plot(xi,[S.lambdify(x,palpha)(i) for i in xi])
        ax.plot(x_samples[0],1.1,'o',ms=10.)
    else:
        g = S.prod(map(f,x_samples[:n]))
        mxval = S.lambdify(a,g)(xi).max()*1.1
        ax.plot(xi,S.lambdify(a,g)(xi))
        ax.plot(x_samples[:n],mxval*ones(n),'o',color='gray',alpha=0.3)
        ax.plot(x_samples[n],mxval,'o',ms=10.)
        ax.axis(ymax=mxval*1.1)
        ax.set_yticklabels([''])
        ax.tick_params(labelsize=6)

# use random  order for first 12 samples
x_samples[:12]= np.random.permutation(x_samples[:12])
fig2,axs = subplots(2,6,sharex=True)
fig2.set_size_inches((10,3))

for n,ax in enumerate(axs.flatten()):
    if n==0:
        ax.plot(xi,[S.lambdify(x,palpha)(i) for i in xi])
        ax.plot(x_samples[0],1.1,'o',ms=10.)
    else:
        g = S.prod(map(f,x_samples[:n]))
        mxval = S.lambdify(a,g)(xi).max()*1.1
        ax.plot(xi,S.lambdify(a,g)(xi))
        ax.plot(x_samples[:n],mxval*ones(n),'o',color='gray',alpha=0.3)
        ax.plot(x_samples[n],mxval,'o',ms=10.)
        ax.axis(ymax=mxval*1.1)
        ax.set_yticklabels([''])
        ax.tick_params(labelsize=6)

fig

import pymc as pm

with pm.Model() as model:
    th = pm.Uniform('th',lower=-pi/2,upper=pi/2)
    alpha= pm.Uniform('a',lower=-3,upper=3)
    xk = pm.Cauchy('x',beta=1,alpha=alpha,observed=x_samples)
    start = pm.find_MAP()
    step = pm.NUTS()
    trace = pm.sample(300, step, start)


pm.traceplot(trace);