Análisis de Referencia

Macros

dots=Dashing[{0.005,0.005}];
ϵϵ=(1/(10^6));app[x_,d_]:=N[((Round[(10^d) x])/(10^d))];
no[x_,μ_,σ_]:=PDF[NormalDistribution[μ,σ],x];
be[x_,α_,β_]:=PDF[BetaDistribution[α,β],x];
ga[x_,α_,β_]:=PDF[GammaDistribution[α,(1/β)],x];
st[x_,θ_,ω_,α_]:=((PDF[StudentTDistribution[α],((x-θ)/ω)])/ω);

Datos Binomiales

Modelo

model[x_,θ_]:=(θ^x)(((1-θ))^(1-x));
model[x,θ]

randata[n_,θ_]:=RandomInteger[BernoulliDistribution[θ],n];
suf[data_]:=Sum[data[[i]],{i,1,Length[data]}]
θθ=0.3;nn=20;
dd=randata[nn,θθ];
rr=suf[dd]; {rr,nn}

Mean[BernoulliDistribution[θ]]
Variance[BernoulliDistribution[θ]]
Mean[BinomialDistribution[n,θ]]
Variance[BinomialDistribution[n,θ]]

FullSimplify[ExpectedValue[(((x- n θ))^2), BinomialDistribution[n,θ],x],{n∈Integers, n>0,0<θ<1}]

La funcion de verosimilitud es
(∏_(i=1)^n)(θ^(Subscript[x, i]))(((1-θ))^(1-(Subscript[x, i])))  =(θ^r)(((1-θ))^(n-r))

like[θ_,r_,n_]:=(θ^r) (((1-θ))^(n-r));
{rr,nn}
Plot[like[θ,rr,nn],{θ,0,1},PlotRange->All]

Solve[D[Log[like[θ,r,n]],θ]==0,θ]

FullSimplify[D[Log[like[θ,r,n]],{θ,2}]/.{θ->(r/n)},{0≤r≤n,n>0, n∈Integers}]

mlest[r_,n_]:=(r/n);
{θθ, N[mlest[rr,nn]]}

Función inicial de referencia

Prior (θ)=(Sqrt[i(θ)])
i(θ)= - E(D(log p(x/θ),(θ,2)))

fisherdata=FullSimplify[-D[Log[model[x,θ]],{θ,2}],{x∈{0,1}, 0<θ<1}]

FullSimplify[ExpectedValue[fisherdata, BernoulliDistribution[θ],x],{0<θ<1}]

FullSimplify[fisherdata/.{x->θ},{0<θ<1}]

i[θ_]=(1/(θ-(θ^2)));
FullSimplify[(Sqrt[i[θ]]),{0<θ<1}]

Con datos binomiales, la distribución  inicial de referencia es
π[θ]=(Sqrt[i(θ)])=(θ^(-1/2))(((1-θ))^(-1/2))

prior[θ_]:=(1/(Sqrt[-(-1+θ) θ]));
prior[θ]

Distribución final de referencia

kernel[θ_,r_,n_]:=like[θ,r,n]prior[θ];
FullSimplify[kernel[θ,r,n],{ n>0, 0≤r≤n,0<θ<1}]

crn=Integrate[kernel[θ,r,n],{θ,0,1},Assumptions->{n>0, 0≤r≤n}]

FullSimplify[((kernel[θ,r,n])/crn), {n>0, 0≤r≤n,0<θ<1}]

post[θ_,r_,n_]:=(((((1-θ))^(-(1/2)+n-r)) (θ^(-(1/2)+r)) Gamma[1+n])/(Gamma[(1/2)+n-r] Gamma[(1/2)+r]));

be[θ,r+1/2,n-r+1/2]

?Beta

Con datos  binomiales, la distribución  posterior  de referencia es
π[θ| r, n]∝ (θ^(r-1/2))(((1-θ))^(n-r-1/2))∝ Be(θ| r+1/2, n-r+1/2)

Plot[post[θ,rr,nn],{θ,0,1},PlotRange->All]

Plot[post[θ,rr,nn],{θ,0,1},Frame->True,Axes->False,PlotRange->All]

nn=100;Plot[post[θ,0,nn],{θ,0,1},Frame->True,Axes->False,PlotRange->{0,10}]

Plot[post[θ,nn,nn],{θ,0,1},Frame->True,Axes->False,PlotRange->{0,10}]

P[a<θ<b | r, n]= (∫_a^b)π[θ| r, n]dθ=(∫_a^b) Be(θ| r+1/2, n-r+1/2)dθ

prob[a_,b_,r_,n_]:=NIntegrate[post[θ,r,n],{θ,a,b}];

f1=Plot[post[θ,30,100],{θ,0,1},PlotRange->All]

prob[0.3,1.0,30,100]

f2=Plot[post[θ,30,100],{θ,0.3,1},PlotRange-> All,Filling->Axis];
Show[f1,f2]

f1=Plot[post[θ,0,1000],{θ,0,0.005}]

prob[0.001,1,0,1000]
f2=Plot[post[θ,0,1000],{θ,0.001,0.005},DisplayFunction->Identity,Filling->Axis];
Show[f1,f2,PlotRange->{0,600}]

prob[0.0001,1.0,0,10000]

prob[0.00005,1.0,0,20000]

Datos Exponenciales

Modelo

PDF[ExponentialDistribution[θ],x]

model[x_,θ_]:=(E^(-x θ)) θ;
model[x,θ]

randata[n_,θ_]:=RandomReal[ExponentialDistribution[θ],n];
suf[data_]:=Sum[data[[i]],{i,1,Length[data]}]
θθ=0.5;nn=10;
dd=randata[nn,θθ]
rr=suf[dd]; {rr,nn}

Mean[ExponentialDistribution[θ]]
Variance[ExponentialDistribution[θ]]

La funcion de verosimilitud es
(∏_(i=1)^n)(e^(-(Subscript[x, i]) θ)) θ  =(θ^n)(e^(- r θ))

like[θ_,r_,n_]:=(θ^n) (E^(-r θ));
Plot[like[θ,rr,nn],{θ,0,2},PlotRange->All]

Solve[D[Log[like[θ,r,n]],θ]==0,θ]

FullSimplify[D[Log[like[θ,r,n]],{θ,2}]/.{θ->(n/r)},{r>0,n>0}]

mlest[r_,n_]:=(n/r);
{θθ, N[mlest[rr,nn]]}

Función inicial de referencia

fisherdata=FullSimplify[-D[Log[model[x,θ]],{θ,2}],{x>0, θ>0}]

i[θ_]=(1/(θ^2));
FullSimplify[(Sqrt[i[θ]]),{θ>0}]

Con  datos exponenciales,la distribución  inicial  de  referencia es
π[θ]=(Sqrt[i(θ)])=(θ^(-1))

prior[θ_]:=(θ^(-1));
prior[θ]

Distribución final de referencia

kernel[θ_,r_,n_]:=like[θ,r,n]prior[θ];
FullSimplify[kernel[θ,r,n],{ n>0, r>0,θ>0}]

crn=Integrate[kernel[θ,r,n],{θ,0,∞},Assumptions->{n>0, r>0}]

FullSimplify[((kernel[θ,r,n])/crn), {n>0,  r>0,θ>0}]

post[θ_,r_,n_]:=(((E^(-r θ)) (r^n) (θ^(-1+n)))/(Gamma[n]));
post[θ,r,n]

ga[θ,n,r]

Con  datos  exponenciales, la distribución  posterior  de referencia es
π[θ| r, n]∝(e^(-r θ))  (θ^(n-1))∝ Ga(θ| n, r)

Plot[post[θ,rr,nn],{θ,0,2},PlotRange->All]

Plot[post[θ,20,100],{θ,0,20},PlotRange->All]

P[a<θ<b | r, n]= (∫_a^b)π[θ| r, n]dθ=(∫_a^b) Ga(θ|n, r)dθ

prob[a_,b_,r_,n_]:=NIntegrate[post[θ,r,n],{θ,a,b}];

f1=Plot[post[θ,rr,nn],{θ,0,2},PlotRange->All]

prob[0.5,∞,rr,nn]
f2=Plot[post[θ,rr,nn],{θ,0.5,2},PlotRange->All,Filling->Axis];
Show[f1,f2]

f1=Plot[post[θ,1,1000],{θ,800,1200},PlotRange->All]

prob[950,1050,1,1000]
f2=Plot[post[θ,1,1000],{θ,950,1050},PlotRange-> All,Filling->0];
Show[f1,f2,PlotRange->All]

Distribución Predictiva de Referencia

Integrate[model[x,θ]post[θ,r,n],{θ,0,∞},Assumptions->{x>0,r>0,n>0}]

pred[x_,r_,n_]:=n (r^n) (((r+x))^(-1-n));
{rr,nn}
f3=Plot[pred[x,rr,nn],{x,0,20},PlotRange->All]

f4=Plot[model[x,θθ],{x,0,20},PlotRange->All,PlotStyle->Red];
Show[f3,f4,PlotRange->All]

Integrate[pred[x,r,n],{x,t,∞},Assumptions->{r>0,t>0,n>0}]

probpred[t_,r_,n_]:=((((r/(r+t))))^n);
probpred[5,rr,nn]//N
f5=Plot[pred[x,rr,nn],{x,5,20},PlotRange->All,Filling->Axis];
Show[f3,f5,PlotRange->All]

Datos normales

Modelo

PDF[NormalDistribution[μ,σ],x]

model[x_,μ_,σ_]=((E^(-((((x-μ))^2)/(2 (σ^2)))))/((Sqrt[2 π]) σ));
model[x,μ,σ]

randata[n_,μ_,σ_]:=RandomReal[NormalDistribution[μ,σ],n];
suf[data_]:={Mean[data],StandardDeviation[data]}
μμ=5;σσ=0.5;nn=10;
dd=randata[nn,μμ,σσ]
mm=suf[dd][[1]]; ss=suf[dd][[2]]; {nn,mm,ss}

Mean[NormalDistribution[μ,σ]]
Variance[NormalDistribution[μ,σ]]

La función de verosimilitud es
(∏_(i=1)^n)N((Subscript[x, i]) |  μ, σ)∝ (σ^(-n))Exp[-(n/(2 (σ^2)))((s^2)+ ((((OverBar[x])-μ))^2))]

like[μ_,σ_,m_,s_,n_]:=(σ^(-n)) (E^(-((n ((s^2)+(((m-μ))^2)))/(2 (σ^2)))));
Plot3D[like[μ,σ,mm,ss,nn],{μ,4,6},{σ,0.1,1.5},PlotPoints->30,PlotRange->All]
ContourPlot[like[μ,σ,mm,ss,nn],{μ,4,6},{σ,0.1,1.5},PlotPoints->100,PlotRange->All]

Solve[{D[Log[like[μ,σ,m,s,n]],μ]==0,D[Log[like[μ,σ,m,s,n]],σ]==0},{μ,σ}]

{{μμ,σσ},{mm,ss}}

Función inicial de referencia

fisherdata=FullSimplify[{
{-D[Log[model[x,μ,σ]],{μ,2}],-D[D[Log[model[x,μ,σ]],μ],σ]},
{-D[D[Log[model[x,μ,σ]],σ],μ],-D[Log[model[x,μ,σ]],{σ,2}]}},{x∈Reals, μ∈Reals, σ>0}];
MatrixForm[fisherdata]

fisher=FullSimplify[ExpectedValue[fisherdata,NormalDistribution[μ,σ],x],{x∈Reals, μ∈Reals, σ>0}];
MatrixForm[fisher]

priormarg[μ_]:=1;
priorcond[σ_,μ_]:=(1/σ);
prior[μ_,σ_]:=priorcond[σ,μ]priormarg[μ];
prior[μ,σ]

Con  datos  normales, la distribución  inicial de referencia para μ 
π[μ, σ]=(σ^(-1))
Analogamente, la distribución  inicial de referencia para σ es tambien
π[μ, σ]=(σ^(-1))

Distribución final de referencia para μ

π[μ| m,s, n]∝ (∫_0^∞)like[μ,σ,m,s,n] prior[μ,σ]dσ

Integrate[like[μ,σ,m,s,n]prior[μ,σ],{σ,0,∞}, Assumptions->{μ∈Reals, n>1, m∈Reals, s>0}]

Posterior proporcional a 
((((s^2)+(((m-μ))^2)))^(-n/2 ))

FullSimplify[st[μ,m,(s/(Sqrt[n-1])),n-1],{n>1, m∈Reals, s>0,μ∈Reals}]

post[μ_,m_,s_,n_]:=st[μ,m,(s/(Sqrt[n-1])),n-1];
{mm,ss,nn}
Plot[post[μ,mm,ss,nn],{μ,4,6}]

nn=200;dd1=randata[nn,300,0.1];
mm=suf[dd1][[1]];ss=suf[dd1][[2]];{nn,mm,ss}
aa=299.95;bb=300.05;
f1=Plot[post[μ,mm,ss,nn],{μ,aa,bb},PlotRange->All]

nn=200;dd2=randata[nn,300,0.1];
mm=suf[dd2][[1]];ss=suf[dd2][[2]];{nn,mm,ss}
f2=Plot[post[μ,mm,ss,nn],{μ,aa,bb},PlotRange->All]

nn=200;dd3=randata[nn,300,0.1];
mm=suf[dd3][[1]];ss=suf[dd3][[2]];{nn,mm,ss}
f3=Plot[post[μ,mm,ss,nn],{μ,aa,bb},PlotRange->All]

Show[f1,f2,f3]

dd=Join[dd1,dd2,dd3];
mm=suf[dd][[1]];ss=suf[dd][[2]];nn=Length[dd];{nn,mm,ss}
ft=Plot[post[μ,mm,ss,nn],{μ,299.95`,aa,bb},PlotRange->All,PlotStyle->Red];
Show[f1,f2,f3,ft]

Distribución final de referencia para σ

π[σ| m,s, n]∝ (∫_(-∞)^∞)like[μ,σ,m,s,n] prior[μ,σ]dμ

Integrate[like[μ,σ,m,s,n]prior[μ,σ],{μ,-∞,∞}, Assumptions->{σ>0, n>1, m∈Reals, s>0}]

cts=Integrate[(((E^(-((n (s^2))/(2 (σ^2))))) (Sqrt[2 π]) (σ^(-n)))/(Sqrt[n])),{σ,0,∞},Assumptions->{n>1,  s>0}]

FullSimplify[((((E^(-((n (s^2))/(2 (σ^2))))) (Sqrt[2 π]) (σ^(-n)))/(Sqrt[n]))/cts),{n>1,  s>0,σ>0}]

postσ[σ_,s_,n_]:=(((2^((3/2)-(n/2))) (E^(-((n (s^2))/(2 (σ^2))))) (((n (s^2)))^((1/2) (-1+n))) (σ^(-n)))/(Gamma[(1/2) (-1+n)]));

μμ=5;σσ=0.5;nn=10;
dd=randata[nn,μμ,σσ];
mm=suf[dd][[1]];ss=suf[dd][[2]];{nn,mm,ss}
Plot[postσ[σ,ss,nn],{σ,0.01`,2.`},PlotRange->All]

μμ=5;σσ=0.5;nn=100;
dd=randata[nn,μμ,σσ];
mm=suf[dd][[1]];ss=suf[dd][[2]];{nn,mm,ss}
Plot[postσ[σ,ss,nn],{σ,0.01`,2.`},PlotRange->All]

Ejemplo numerico

http://en.wikipedia.org/wiki/Earth%27s_gravity
Madrid 9.800 m/s2.  Paris 9.809 m/s2

dd={9.8,9.796,9.81,9.802,9.816}
mm=suf[dd][[1]];ss=suf[dd][[2]];nn=Length[dd];
{nn,mm,ss}
Plot[post[μ,mm,ss,nn],{μ,9.76,9.84},Frame->True,Axes->False,PlotRange->All]
Plot[postσ[σ,ss,nn],{σ,0,0.05`},Frame->True,Axes->False,PlotRange->All]

NIntegrate[post[μ,mm,ss,nn],{μ,9.81,∞}]
NIntegrate[post[μ,mm,ss,nn],{μ,9.82,∞}]
NIntegrate[postσ[σ,ss,nn],{σ,0.01,∞}]
NIntegrate[postσ[σ,ss,nn],{σ,0.02,∞}]