O Limite Mínimo da Variância

Cálculo teórico para uma f.d.p. de Cauchy

f=PDF[CauchyDistribution[θ,1],x]

Integrate[f,{x,-Infinity,Infinity},Assumptions->θ∈Reals] (* usando a propriedade da normalização para verificar que achamos como integrar essa função *)

LMV é o inverso e oposto do valor esperado da segunda derivada - começamos derivando

d2f=D[Log[f],{θ,2}]

d2f=FullSimplify[d2f]

A função verossimilhança é o produto de n funções de Cauchy e a derivada do LogL dá uma soma de n termos. Assim, terminamos com a soma de n termos identicos...

lmv=-1/(n Integrate[f d2f,{x,-Infinity,Infinity},Assumptions->θ∈Reals])

Comparação com dois estimadores diferentes...

Determinação da variância da mediana como estimativa de θ

simulação...

u=RandomVariate[CauchyDistribution[0,1],10]

Median[u]

repetir a simulação muitas vezes (1000) e determinar a variancia da estimativa, para comparar com o LMV

q=Table[Median[RandomVariate[CauchyDistribution[0,1],10]],{1000}];

Variance[q]

Repetir o cálculo para vários valores de n, pequenos

varMediana=Table[ {n,Variance[Table[Median[RandomVariate[CauchyDistribution[0,1],n]],{10000}]]},{n,3,30}]

Show[{ListPlot[varMediana,PlotRange->All,AxesLabel->{"n","variance"},BaseStyle->{FontSize->14}],Plot[2/n,{n,3,30}]}]

Show[{ListPlot[varMediana,PlotRange->{0,0.5},AxesLabel->{"n","variance"},BaseStyle->{FontSize->14}],Plot[2/n,{n,3,30}]}]

Para n grande, há outra forma de estimar mais precisa que a mediana

Procedimento: 1) O conjunto de dados é ordenado - como estamos simulando, encadeamos as duas operações

u=Sort[RandomVariate[CauchyDistribution[0,1],200]];

e 2) calcula-se a média dos 24% dos dados em torno da mediana; Floor e Ceiling arredondam para Integer

Mean[u[[Floor[0.38 200];;Ceiling[0.62 200]]] ]

Dá para fazer em uma instrução só, embora fique meio ilegível

Mean[Sort[RandomVariate[CauchyDistribution[0,1],200]][[Floor[0.38 200];;Ceiling[0.62 200]]] ]

repete a simulaçao 100 vezes e calcula a variancia da estimativaa

Variance[Table[Mean[Sort[RandomVariate[CauchyDistribution[0,1],200]][[Floor[0.38 200];;Ceiling[0.62 200]]] ],{1000}]]

fazemos agora para diversos valores do número de dados n

varEstim24=Table[ {n,Variance[Table[Mean[Sort[RandomVariate[CauchyDistribution[0,1],n]][[Floor[0.38 n];;Ceiling[0.62 n]]] ],{1000}]]},{n,100,2000,100}]

fazemos o mesmo para a mediana, para comparar

varMedianaNGrande=Table[ {n,Variance[Table[Median[RandomVariate[CauchyDistribution[0,1],n]],{1000}]]},{n,100,2000,100}]

ListPlot[{varMedianaNGrande,varEstim24},PlotRange->All,AxesLabel->{"n","variance"},BaseStyle->{FontSize->14}]

comparacao=ListPlot[{varMedianaNGrande,varEstim24},PlotRange->{0,0.004},PlotLegends->{"mediana","media central"},AxesLabel->{"n","variance"},BaseStyle->{FontSize->14}]

Para que os pontos simulados variem mais suavemente, basta repetir os calculos com tiragens de 10000 simulações...

Show[{comparacao,Plot[2/n,{n,200,2000}]}]

Maximum Likelihood estimate of Cauchy localization parameter

Cauchy distribution location parameter with unity width

preparing / testing

PDF[CauchyDistribution[x0,γ],x]

logLCauchy[data_][a_]:=-Total[Log[1+(data-a)^2] ]

w=RandomVariate[CauchyDistribution[10,1],10]

logLc=logLCauchy[w]

Plot[logLc[a],{a,5,15},AxesLabel->{"a","LogL"},BaseStyle->{FontSize->14}]

estimating location parameter from the maximum

ahat[data_]:=Module[{u,a},u=data-a;
a/.FindRoot[Total[data/(1+u^2)]==a Total[1/(1+u^2)],{a,Median[data]}]
]

ahat[w]

demonstration

Manipulate[randomize; Module[{u,logL,σ0,min,max},u=RandomVariate[CauchyDistribution[a0,1],n] ;
logL=logLCauchy[u];
σ0=Sqrt[2/n];
min=Max[0,a0-5σ0];
max = a0+5σ0;
Plot[{ logL[a],-(a-a0)^2/(2σ0^2)+logL[a0]},{a,min,max},BaseStyle->{FontSize->14},PlotStyle->{Blue,Red},AxesLabel->{a,"logL"},PlotLabels->{"exp","limit"},Epilog->{Dashed,Red,Line[{{a0,logL[max]},{a0,logL[a0]}}]}]
],{{a0,3},1,10},{{n,10},1,100,1},Row[{Button["Sort again",{randomize=RandomReal[]},Appearance->"Palette",ImageSize->{240,28}]}],ControlPlacement->Bottom]

Comparison between the efficiencies of the Median, trimmed median and Maximum Likelihood

varMaxL=Table[ {n,Variance[Table[ahat[RandomVariate[CauchyDistribution[0,1],n]],{100}]]},{n,15,30}]

Show[{ListPlot[{varMediana,varMaxL},PlotRange->{0,0.5},AxesLabel->{"n","variance"},BaseStyle->{FontSize->14}],Plot[2/n,{n,3,30}]}]

Example with n=100

varaHat100ML=Variance[Table[ahat[RandomVariate[CauchyDistribution[0,1],100]],{1000}]]

2/100//N

varaHatTrimMedian=Variance[Table[Mean[Sort[RandomVariate[CauchyDistribution[0,1],100]][[Floor[0.38 100];;Ceiling[0.62 100]]] ],{1000}]]