%matplotlib inline
from __future__ import  division
from IPython.html.widgets import interact,fixed
from scipy import stats
import sympy
from sympy import S, init_printing, Sum, summation,expand
from sympy import stats as st
from sympy.abc import k
import re

sympy.plot( k*(1-k),(k,0,1),ylabel='$n^2 \mathbb{V}$',xlabel='p',fontsize=28)

b=stats.bernoulli(p=.8)
samples=b.rvs(100)
print var(samples)
print mean(samples)

def slide_figure(n=100,m=1000,p=0.8):
    fig,ax=subplots()
    ax.axis(ymax=25)
    b=stats.bernoulli(p=p)
    v=iter(b.rvs,2)
    samples=b.rvs(n*1000)
    tmp=(samples.reshape(n,-1).mean(axis=0))
    ax.hist(tmp,normed=1)
    ax1=ax.twinx()
    ax1.axis(ymax=25)
    ax1.plot(linspace(0,1,200),stats.norm(mean(tmp),std(tmp)).pdf(linspace(0.0,1,200)),lw=3,color='r')

interact(slide_figure, n=(100,500),p=(0.01,1,.05),m=fixed(500));

from sympy.abc import p,n,k

sympy.plot(st.density(st.Beta('p',3,3))(p),(p,0,1) ,xlabel='p')

obj=sympy.expand_log(sympy.log(p**k*(1-p)**(n-k) * st.density(st.Beta('p',6,6))(p)))

sol=sympy.solve(sympy.simplify(sympy.diff(obj,p)),p)[0]
print sol

n=5
def show_bias(n=30):
    sympy.plot(p,(5+n*p)/(n+10),(p,0,1),aspect_ratio=1,title='more samples reduce bias')
    
interact(show_bias,n=(10,500,10));

sum(sympy.var('x1:10'))
expr=((5+(sum(sympy.var('x1:10'))))/(n+10))

def apply_exp(expr):
    tmp=re.sub('x[\d]+\*\*2','p',str(expand(expr)))
    tmp=re.sub('x[\d]+\*x[\d]+','p**2',tmp)
    tmp=re.sub('x[\d]+','p',tmp)
    return sympy.sympify(tmp)

ex2 = apply_exp(expr**2)
print ex2

tmp=sympy.simplify(ex2 - (apply_exp(expr))**2 )
sympy.plot(tmp,p*(1-p)/10,(p,0,1))

def generate_expr(num_samples=10,alpha=6):
    n = sympy.symbols('n')
    obj=sympy.expand_log(sympy.log(p**k*(1-p)**(n-k) * st.density(st.Beta('p',alpha,alpha))(p)))
    sol=sympy.solve(sympy.simplify(sympy.diff(obj,p)),p)[0]
    expr=sol.replace(k,(sum(sympy.var('x1:%d'%(num_samples)))))
    expr=expr.subs(n,num_samples)
    ex2 = apply_exp(expr**2)
    ex =  apply_exp(expr)
    return (ex,sympy.simplify(ex2-ex**2))

num_samples=10
X_bias,X_v = generate_expr(num_samples,alpha=2)
p1=sympy.plot(X_v,(p,0,1),show=False,line_color='b',ylim=(0,.03),xlabel='p')
X_bias,X_v = generate_expr(num_samples,alpha=6)
p2=sympy.plot(X_v,(p,0,1),show=False,line_color='r',xlabel='p')
p3=sympy.plot(p*(1-p)/num_samples,(p,0,1),show=False,line_color='g',xlabel='p')
p1.append(p2[0])
p1.append(p3[0])
p1.show()

p1=sympy.plot(n*(1-p)*p/(n+10)**2,(p,0,1),show=False,line_color='b',ylim=(0,.05),xlabel='p',ylabel='variance')
p2=sympy.plot((1-p)*p/n,(p,0,1),show=False,line_color='r',ylim=(0,.05),xlabel='p')
p1.append(p2[0])
p1.show()

def show_variance(n=5):
    p1=sympy.plot(n*(1-p)*p/(n+10)**2,(p,0,1),show=False,line_color='b',ylim=(0,.05),xlabel='p',ylabel='variance')
    p2=sympy.plot((1-p)*p/n,(p,0,1),show=False,line_color='r',ylim=(0,.05),xlabel='p')
    p1.append(p2[0])
    p1.show()    
interact(show_variance,n=(10,120,2));


pv=.60
nsamp=30
fig,ax=subplots()
rv = stats.bernoulli(pv)
map_est=(rv.rvs((nsamp,1000)).sum(axis=0)+5)/(nsamp+10);
ml_est=(rv.rvs((nsamp,1000)).sum(axis=0))/(nsamp);
_,bins,_=ax.hist(map_est,bins=20,alpha=.3,normed=True,label='MAP');
ax.hist(ml_est,bins=20,alpha=.3,normed=True,label='ML');
ax.vlines(map_est.mean(),0,12,lw=3,linestyle=':',color='b')
ax.vlines(pv,0,12,lw=3,color='r',linestyle=':')
ax.vlines(ml_est.mean(),0,12,lw=3,color='g',linestyle=':')
ax.legend()