Método da máxima verossimilhança no ajuste de uma reta, dados \[Distributed] Poisson

Exemplo: espectro de energia em que o número de contagens y(x) = a + b x, onde x é um número inteiro, a ~1 e b~0

Este é o conjunto de dados, tirado de um experimento e relocado por questões didáticas

dados={{11,0},{12,3},{13,4},{14,5},{15,6},{16,5},{17,0},{18,3},{19,2},{20,0},{21,1},{22,1},{23,2},{24,2},{25,1},{26,3}}

ListPlot[dados]

Needs["ErrorBarPlots`"]

ErrorListPlotPoisson[data_]:=Module[{u,t},u=N[Sqrt[data[[All,2]] ] ];
t=Partition[ Flatten[Transpose[{dados,u/.x_/;x<1->1}] ],3];
ErrorListPlot[t ,PlotRange->All]]

ErrorListPlotPoisson[dados]

solução gráfica

logL[data_,a_,b_]:=-Length[data]a-Total[data[[All,1]] ]b+data[[All,2]].Log[a+b×data[[All,1]] ]

ContourPlot[Evaluate[logL[dados,a,b] ],{a,2,8},{b,-0.2,0}]

gradeLarga=Table[{a,b,logL[dados,a,b]} ,{a,4,6,0.4},{b,-0.2,-0.05,0.03}]

gradeLarga//TableForm

Max[ Select[ Flatten[gradeLarga,1][[All,3]],# ∈ Reals&] ]

gradeMedia=Table[{a,b,logL[dados,a,b]} ,{a,4,4.8,0.2},{b,-0.14,-0.08,0.01}]

TableForm[gradeMedia]

gradeFina=Table[{a,b,logL[dados,a,b]} ,{a,4.4,4.8,0.05},{b,-0.13,-0.11,0.0025}]

TableForm[gradeFina]

logL = máximo - n^2/2 corresponde a n desvios-padrões quando há só uma variável, aproximadamente

ContourPlot[Evaluate[logL[dados,a,b] ],{a,0,10},{b,-0.3,0.1}, Contours->{-4.5,-6,-8.5}]

Solução analítica

eqA=Length[dados]==Total[dados[[All,2]]/(ap+bp dados[[All,1]] ) ]

eqB=Total[ dados[[All,1]] ]==Total[dados[[All,2]]dados[[All,1]]/(ap+bp dados[[All,1]] ) ]

solucoes=Solve[{eqA,eqB},{ap,bp}]//N  (* Mathematica tem uma representação exata para estas soluções, porque não envolve funções transcendentes em ap e bp, mas são complicadas demais para serem úteis, por isso ficamos com as soluções numéricas *)

é preciso selecionar as raízes que dão logL real

logLcalc=Table[logL[dados,ap,bp]/.solucoes[[i]],{i,Length[solucoes]}]

TableForm[Transpose[{solucoes,logLcalc}],TableHeadings->{None,{"ponto","logL"}}]

dessa tabela, verifica-se que solucoes[[7]] é a raiz - vamos guardar numa variável para uso futuro

solucao=solucoes[[7]]

ErrorListPlotPoisson[dados]

Show[{ErrorListPlotPoisson[dados],Plot[ap+bp x/.solucao,{x,10,26}]}]

agora, as barras de incerteza podem ser calculadas com uma aproximação melhor

Show[{ErrorListPlot[Partition[ Flatten[Transpose[{dados,Sqrt[ap+bp dados[[All,1]]/.solucao]}]],3] ],Plot[ap+bp x/.solucao,{x,10,26}]}]

residuo=dados[[All,2]]-(ap+bp dados[[All,1]]/.solucao)

residuo/Sqrt[ap+bp dados[[All,1]]/.solucao]

Sempre graficar os residuos normalizados em funçao das variáveis dependente e independente

ListPlot[ Transpose[{dados[[All,1]],residuo/Sqrt[ap+bp dados[[All,1]]/.solucao] }]]

7/16//N

usamos o valor ajustado para a coordenada independente

ListPlot[ Transpose[{(ap+bp dados[[All,1]]/.solucao),residuo/Sqrt[ap+bp dados[[All,1]]/.solucao] }]]

variância a partir da matriz das derivadas da função log(verossimilhança) - exige o valor da estimativa!

ycalc=ap+bp dados[[All,1]]/.solucao

matrixM={{Total[dados[[All,2]]/ycalc^2],Total[dados[[All,1]] dados[[All,2]]/ycalc^2]},{Total[dados[[All,1]] dados[[All,2]]/ycalc^2],Total[dados[[All,1]]^2 dados[[All,2]]/ycalc^2]} }

varAB=Inverse[matrixM]

desvioPadrao=Round[Sqrt[Diagonal[varAB] ],0.01]

solucao

Resultado: a = 4.6(16), b = -0.12(8)

determinação da tendenciosidade por simulação

simul[a_,b_,x_]:=Table[RandomVariate[PoissonDistribution[a+b x[[i]] ] ],{i,Length[x]}]

dadosSimulados=Transpose[{Range[11,26],simul[4,-0.1,Range[11,26]]}]

o método analítico é mais indicado para repetir muitas vezes o cálculo

equações de máximo da verossimilhança

eqAp[dados_]:=Length[dados]==Total[dados[[All,2]]/(ap+bp dados[[All,1]] ) ]

eqBp[dados_]:=Total[ dados[[All,1]] ]==Total[dados[[All,2]]dados[[All,1]]/(ap+bp dados[[All,1]] ) ]

resolve o sistema

solucoes=Solve[{eqAp[dadosSimulados],eqBp[dadosSimulados]},{ap,bp}]//N (* esse //N transforma em númerico os resultados analíticos encontrados pelo Mathematica, uma vez que as equações em ap e bp não envolvem nenhuma função transcedente - Mathematica tem uma solução exata para esse caso, mas que é complicada demais para ser útil *)

seleção das raízes em duas etapas: primeiro isola as reais,

logLcalc=Table[logL[dadosSimulados,ap,bp]/.solucoes[[i]],{i,Length[solucoes]}]

Select[logLcalc,# ∈ Reals &]

e depois separa a que dá o máximo de logL

Position[logLcalc,Max[Select[logLcalc,# ∈ Reals &] ] ]

solucoes[[Flatten[Position[logLcalc,Max[Select[logLcalc,# ∈ Reals &] ] ] ] ]]

compacta toda a solução em uma função

solve[dados_]:=Module[{solucoes,logLcalc},
solucoes=Solve[{eqAp[dados],eqBp[dados]},{ap,bp}]//N;
logLcalc=Table[logL[dados,ap,bp]/.solucoes[[i]],{i,Length[solucoes]}];
solucoes[[Flatten[Position[logLcalc,Max[Select[logLcalc,# ∈ Reals &] ] ] ] ]][[1]]
]

solve[dadosSimulados]

prepara a simulação

extrai as coordenadas x de interesse a partir dos dados

xexp=dados[[All,1]]

é mais prático criar uma função que gera a tabela {x,y}, com y simulado

simulDados[a_,b_,x_]:=Table[{x[[i]],RandomVariate[PoissonDistribution[a+b x[[i]] ] ]},{i,Length[x]}]

simulDados[4,-0.1,xexp]

simulando...

solucao

ap/.solucao

começar com 100 replicas e ir aumentando aos poucos, se há tempo!

replicas=Table[solve[simulDados[ap/.solucao,bp/.solucao,xexp ] ],{100}]

Mean[{ap,bp}/.replicas]

StandardDeviation[{ap,bp}/.replicas]/Sqrt[100]

não dá para ter certeza que é tendencioso, valores conferem em menos que 2 desvios padroes... tentar 1000 repeticoes

replicas=Table[solve[simulDados[4.6,-0.115,xexp ] ],{1000}];

Mean[{ap,bp}/.replicas]

desvio=StandardDeviation[{ap,bp}/.replicas]

desvio/Sqrt[1000]

diferenca=({ap,bp}/.solucao)-Mean[{ap,bp}/.replicas]

Conclusão: parece que os parâmetros estão ligeiramente superestimados, mas a tendenciosidade é muito menor que os desvios-padrões do resultado, de modo que não vale a pena corrigir.

Tendenciosidade muito menor que incerteza, resultado

a = 4,6(17)  e  b = -0,12(8)

Comparando com MMQ, já que a função é linear

ESSENCIAL evitar pesos nulos no MMQ

pesos=1/(dados[[All,2]]/.x_/;x<1->1)

mmqIncompleto=LinearModelFit[dados,{x},x,Weights->pesos,VarianceEstimatorFunction->(1&)]

mmqIncompleto["ParameterTable"]

abI=mmqIncompleto["ParameterTableEntries"][[All,1]]

Show[{ErrorListPlot[Partition[ Flatten[Transpose[{dados,Sqrt[abI[[1]]+abI[[2]] dados[[All,1]]/.solucao]}]],3] ],Plot[abI[[1]]+abI[[2]] x,{x,9,26}]}]

Esse resultado é muito ruim, mas ele vem da subestimação da variância por conta da assimetria da f.p. de Poisson. Por isso, fazemos outro ciclo...

pesos = 1/(abI[[1]]+abI[[2]] dados[[All,1]])

mmqI1=LinearModelFit[dados,{x},x,Weights->pesos,VarianceEstimatorFunction->(1&)]

mudou muito - iterar outra vez, com novos pesos

abI=mmqI1["ParameterTableEntries"][[All,1]]

pesos = 1/(abI[[1]]+abI[[2]] dados[[All,1]])

mmqI2=LinearModelFit[dados,{x},x,Weights->pesos,VarianceEstimatorFunction->(1&)]

mais uma iteração ...

abI=mmqI2["ParameterTableEntries"][[All,1]]

mmqI3=LinearModelFit[dados,{x},x,Weights->1/(abI[[1]]+abI[[2]] dados[[All,1]]),VarianceEstimatorFunction->(1&)]

uma última iteração, por desencargo de consciência...

abI=mmqI3["ParameterTableEntries"][[All,1]]

mmq=LinearModelFit[dados,{x},x,Weights->1/(abI[[1]]+abI[[2]] dados[[All,1]]),VarianceEstimatorFunction->(1&)]

mmq["ParameterTable"]

compare com resultado da máxima verossimilhança: a = 4,6(17)  e  b = -0,12(9)

Assim, NESTE caso em que a função é linear nos parâmetros, o mmq se sai muito bem, DESDE QUE se refaça os cálculos com pesos deduzidos dos valores ajustados, em um processo iterativo, até a convergência. Usar Sqrt[ n experimental ] para o desvio-padrão de uma grandeza com f.p. Poisson é uma aproximação inadequada quando n < 10 e muito ruim quando n ~ 1.