Integración por Montecarlo

Parte 1 : aproximación al seno mediante una uniforme

resultado[f_,a_,b_]:=Integrate[f[x],{x,a,b},Assumptions->{x∈Reals,f[x]≥0,a<b}];

ff[x_]:=Sin[x];

aa=0;bb=π;

Verdadero valor

rr=resultado[ff,aa,bb]

Gráficas

f1=Plot[ff[x],{x,aa,bb},PlotRange->All]

pp[x_,a_,b_]:=If[a<x<b,(1/(b-a)),0];

f2=Plot[pp[x,aa,bb],{x,aa,bb},PlotStyle->Red,PlotRange->All];

Show[f1,f2]

Generamos números al azar de la uniforme

rmd[n_,a_,b_]:=Table[RandomReal[{a,b}],{i,1,n}];
rmd[10,0,20]

Aproximación numérica al resultado real

intaprox[f_,a_,b_,p_,n_]:=Module[{y},y=rmd[n,a,b];
(1/n)Sum[((f[y[[i]]])/(p[y[[i]],a,b])),{i,1,n}]];

nn=10;intaprox[ff,aa,bb,pp,nn]

nn=10000;intaprox[ff,aa,bb,pp,nn]

nn1=5;nn2=100000;dd=1000;
resultados=Table[{n,intaprox[ff,aa,bb,pp,n]},{n,nn1,nn2,dd}];
g0=Plot[rr,{n,nn1,nn2},PlotStyle->Green];
g1=ListPlot[resultados,Joined->True];

hh=0.02;
Show[g0,g1,PlotRange->{rr-hh,rr+hh}]
npoints=Dimensions[resultados][[1]];
last=resultados[[npoints]];
{last,rr-last[[2]]}

nn1=100;nn2=10000;dd=200;hh=0.2;
resultados=Table[{n,intaprox[ff,aa,bb,pp,n]},{n,nn1,nn2,dd}];
g0=Plot[rr,{n,nn1,nn2},PlotStyle->Green];
g1=ListPlot[resultados,Joined->True];
Show[g0,g1,PlotRange->{rr-hh,rr+hh}]
npoints=Dimensions[resultados][[1]];
last=resultados[[npoints]];
{last,rr-last[[2]]}

nn1=100;nn2=100000;dd=10000;hh=0.1;
resultados=Table[{n,intaprox[ff,aa,bb,pp,n]},{n,nn1,nn2,dd}];
g0=Plot[rr,{n,nn1,nn2},PlotStyle->Green];
g1=ListPlot[resultados,Joined->True];
Show[g0,g1,PlotRange->{rr-hh,rr+hh}]
npoints=Dimensions[resultados][[1]];
last=resultados[[npoints]];
{last,rr-last[[2]]}

Modelizacion

nn=100;mm=20;
resultados=Table[intaprox[ff,aa,bb,pp,nn],{i,1,mm}];
media=Mean[resultados]; desv=StandardDeviation[resultados];
{nn,mm,media,desv}

nn=100;mm=10000;
resultados=Table[intaprox[ff,aa,bb,pp,nn],{i,1,mm}];
media=Mean[resultados]; desv=StandardDeviation[resultados];
{nn,mm,media,desv}

No[x_,μ_,σ_]:=PDF[NormalDistribution[μ,σ],x];

Needs["Histograms`"]
fig1=Histogram[resultados,HistogramScale->1];
fig2=Plot[No[x,media,desv],{x,media-5desv, media+5desv}];
Show[fig1,fig2]

nn=10;mm=10000;
resultados=Table[intaprox[ff,aa,bb,pp,nn],{i,1,mm}];
media=Mean[resultados]; desv=StandardDeviation[resultados];
{nn,mm,media,desv}
fig1=Histogram[resultados,HistogramScale->1];
fig2=Plot[No[x,media,desv],{x,media-5desv, media+5desv}];
Show[fig1,fig2]

Hacer inferencia

St[x_,μ_,σ_,α_]:=(1/σ)PDF[StudentTDistribution[α],((x-μ)/σ)];

Plot[St[x,6,1,4],{x,0,12},PlotRange->All]

posterior[μ_,med_,sd_,monts_]:=St[μ,med,(sd/(Sqrt[monts-1])),monts-1];

posterior2[μ_,med_,sd_,monts_]:=If[monts<1000,posterior[μ,med,sd,monts],No[μ,med,(sd/(Sqrt[monts]))]];

nn=100;mm=1000;
Timing[data=Table[intaprox[ff,aa,bb,pp,nn],{j,1,mm}];]
Mean[data]

Integrales por montecarlo

nn=100;mm=20;
Timing[data=Table[intaprox[ff,aa,bb,pp,nn],{j,1,mm}];]
mme=Mean[data];
ssd=StandardDeviation[data];
{mm,mme,ssd}
r0=Plot[posterior[μ,mme,ssd,mm],{μ,mme-((6ssd)/(Sqrt[mm])),mme+((6ssd)/(Sqrt[mm]))},PlotRange->All,PlotStyle->Blue]

nn=1000;mm=2000;
Timing[data=Table[intaprox[ff,aa,bb,pp,nn],{j,1,mm}];]
mme=Mean[data];
ssd=StandardDeviation[data];
{mm,mme,ssd}
r1=Plot[posterior[μ,mme,ssd,mm],{μ,mme-((6ssd)/(Sqrt[mm])),mme+((6ssd)/(Sqrt[mm]))},PlotRange->All, PlotStyle->Green]

nn=100000;mm=20;
Timing[data=Table[intaprox[ff,aa,bb,pp,nn],{j,1,mm}];]
mme=Mean[data];
ssd=StandardDeviation[data];
{mm,mme,ssd}
r2=Plot[posterior[μ,mme,ssd,mm],{μ,mme-((6ssd)/(Sqrt[mm])),mme+((6ssd)/(Sqrt[mm]))},PlotStyle->Red];
Show[r1,r2]

nn=1000000;mm=2;
Timing[data=Table[intaprox[ff,aa,bb,pp,nn],{j,1,mm}];]
mme=Mean[data];
ssd=StandardDeviation[data];
{mm,mme,ssd}
r3=Plot[posterior[μ,mme,ssd,mm],{μ,mme-((6ssd)/(Sqrt[mm])),mme+((6ssd)/(Sqrt[mm]))},PlotStyle->Blue];
Show[r1,r2,r3]

nn=10;mm=200000;
Timing[data=Table[intaprox[ff,aa,bb,pp,nn],{j,1,mm}];]
mme=Mean[data];
ssd=StandardDeviation[data];
{mm,mme,ssd}
r4=Plot[posterior2[μ,mme,ssd,mm],{μ,mme-((6ssd)/(Sqrt[mm])),mme+((6ssd)/(Sqrt[mm]))},PlotStyle->Black];
Show[r1,r2,r3,r4]

nn=2000;mm=5000;
Timing[data=Table[intaprox[ff,aa,bb,pp,nn],{j,1,mm}];]
mme=Mean[data];
ssd=StandardDeviation[data];
{mm,mme,ssd}
r5=Plot[posterior2[μ,mme,ssd,mm],{μ,mme-((6ssd)/(Sqrt[mm])),mme+((6ssd)/(Sqrt[mm]))},PlotStyle->Yellow];
Show[r1,r2,r3,r4,r5]