Método da máxima verossimilhança - 18/11/2014, Vito R. Vanin - Instituto de Física, USP

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,AxesLabel->{c,"contagens/canal"}]

Needs["ErrorBarPlots`"]

juntar uma coluna a uma matriz não é direto

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

ErrorListPlotPoisson[dados]

solução gráfica

logL = máximo - n^2/2 correspondem a n desvios-padrões, 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

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 p=ficamos com as soluções numéricas *)

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

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

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

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

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

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!

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]} }

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

determinação da tendenciosidade por simulação

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

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

resolve o sistema

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

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

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]]
]

prepara a simulação

extrai as coordenadas x de interesse a partir dos dados

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

simulando...

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

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

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

Show[{ErrorListPlot[Join[dados,Transpose[{Sqrt[abI[[1]]+abI[[2]] dados[[All,1]]/.solucao]}],2] ],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...

mudou muito - iterar outra vez, com novos pesos

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.

MMQ para funções não-lineares em apenas 1 ou 2 parâmetros e lineares nos demais

não-linear em um único parâmetro

dados experimentais

seno={{0,490,210},{33,1400,400},{67,-850,350},{100,-130,220}}

graficoDados=ErrorListPlot[seno,PlotRange->All]

Ajuste dos parâmetros

Show[{graficoDados,Plot[fitSemCuidado[t],{t,0,100}]}]

Show[{graficoDados,Plot[fitSemPesos[t],{t,0,100}]}]

Note que a principal mudança por usar os pesos corretos está no desvio-padrão de ω

fitCerto["ParameterTable"]

Verificando a tendenciosidade por simulação

simulacaoMedicao[dados_,a_,ω_]:=Module[{x,σ,y},x=dados[[All,1]];
σ=dados[[All,3]];
Table[{x[[i]],RandomVariate[NormalDistribution[a Sin[ω x[[i]] ],σ[[i]] ] ]},{i,Length[dados]} ] ]

Criar uma função que ajusta e devolve só os parâmetros ajustados

Atenção, não confundir o parâmetro da função simulaFit - ω0 - com a variável ω da função α Sin[ω t] que entra em NonlinearModelFit

simulaFit[dados_,a0_,ω0_]:=NonlinearModelFit[simulacaoMedicao[dados,a0,ω0],{α Sin[ω t]},{{α,1000},{ω,0.06}},{t},Weights->1/dados[[All,3]]^2,VarianceEstimatorFunction->(1&)]["ParameterTableEntries"][[All,1]]

simulaFit[seno, 1268.,0.0615]

às vezes, o algoritmo de NonlinearModelFit pode não encontrar o mínimo com a precisão estabelecida (quando não se estabelece um AccuracyGoal, a rotina usa AccuracyGoal=1/2 PrecisionGoal que, por default, é MachinePrecision=16 dígitos. Assim, ele está tentando encontrar os valores com 8 dígitos de precisão, que pode ser muito. Nesse caso, pode-se mudar simulaFit para incluir na chamada de NolinearModelFit a opção AccuracyGoal->5, por exemplo.

Neste caso, a tendenciosidade é pequena, muito menor que o desvio-padrão do resultado, e não precisa ser corrigida - há algo errado com os números do texto! Já o desvio-padrão estimado merece ser modificado:

Conclusão: A= = 1270(300), ω=0,0615(20)

Frequentemente é melhor desenhar a curva de soma dos resíduos em funçaõ do parâmetro não-linear

Não linear em dois dos parâmetros

qChapeu[dado_,ω_]:=Module[{u,aChapeu},u=Sin[ω dado[[All,1]] ];
aChapeu=Total[ dado[[All,2]] u/dado[[All,3]]^2 ]/Total[u^2/dado[[All,3]]^2 ];
Total[ (dado[[All,2]]-aChapeu u)^2/dado[[All,3]]^2 ]
]

Plot[qChapeu[seno,ω],{ω,0.001,0.2},AxesLabel->{ω,OverHat[Q]}]

quando não se observa o fenômeno em toda sua faixa de variação, é muito pior

simulSeno[a_,σa_,fperiodo_,n_,periodo_]:=Module[{t,u},t= fperiodo periodo Range[n]/n;
u=a Sin[2Pi t /periodo ];
Table[{t[[i]],RandomVariate[ NormalDistribution[u[[i]],σa  ] ],σa },{i,n}]
]

facil=simulSeno[1000,100,1,4,100]

Plot[qChapeu[facil,ω],{ω,0.001,4Pi/100},AxesLabel->{ω,OverHat[Q]}]

mude a incerteza e a fração do período e veja como o mínimo perde definição

Plot[Evaluate[ qChapeu[simulSeno[1000,30,0.1,4,100],ω] ],{ω,0.001,8Pi/100},AxesLabel->{ω,OverHat[Q]}]

Ajuste dos parâmetros a dados com boa precisão

forma funcional

f[t_,a_,b_,w1_,w2_]:=a Sin[w1 t]+b Cos[w2 t]

dados experimentais

p2={{0,170,10},{5,188,11},{10,213,12},{16,235,12},{21,209,12},{26,199,11},{31,200,10},{37,124,8},{42,91,6},{47,37,5},{52,-29,5},{58,-78,6},{63,-125,8},{68,-190,10},{73,-228,12},{79,-261,14},{84,-287,14},{89,-282,15},{94,-274,14},{99,-212,13},{105,-211,11},{110,-154,9},{115,-99,7},{120,-58,5},{126,5,4}}

ErrorListPlot[p2,PlotRange->All]

ajuste e resultados

fit=NonlinearModelFit[p2[[All,{1,2}]],f[t,a1,a2,w1,w2],{{a1,100},{a2,180},{w1,0.05},{w2,0.03}},t,Weights->1/p2[[All,3]]^2,VarianceEstimatorFunction->(1&)]

fit["ParameterTable"]

chiP2=Total[(fit["FitResiduals"]/p2[[All,3]])^2 ]

estimativas precisam ser precisas

NonlinearModelFit[p2[[All,{1,2}]],f[t,a1,w1,a2,w2],{{a1,100},{w1,180},{a2,0.03},{w2,0.08}},t,Weights->1/p2[[All,3]]^2,VarianceEstimatorFunction->(1&)]

gráficos para avaliação da qualidade do ajuste

Show[{ErrorListPlot[p2,PlotRange->All],Plot[fit[t],{t,0,130}]}]

é muito comum graficar-se o resíduo normalizado pelos desvios-padrões dos dados, mas é melhor usar o desvio-padrão do resíduo calculado com a matriz H, como no próximo gráfico, embora neste caso faça pouca diferença porque temos um bom número de dados, com precisões semelhantes

ListPlot[Transpose[{p2[[All,1]],fit["FitResiduals"]/p2[[All,3]]}] ]

ListPlot[Transpose[{p2[[All,1]],fit["StandardizedResiduals"]}] ]

ListPlot[Transpose[{fit["Response"],fit["StandardizedResiduals"]}] ]

Simulando uma situação difícil, com a mesma função do ajuste acima

simulação

dados0= Table[{i Pi/3,f[i Pi/3,100, 180.,0.055,0.035]},{i,0,100,19}]

dp=Round[ Table[Sqrt[(0.20 dados0[[i,2]])^2 + 200],{i,Length[dados0]}] ]

simulacao=Round[ Table[{dados0[[i,1]],RandomVariate[NormalDistribution[dados0[[i,2]],dp[[i]] ] ],dp[[i]]},{i,Length[dados0]}] ]

ErrorListPlot[simulacao,PlotRange->All]

ajuste e resultados

fitSimulacao=NonlinearModelFit[simulacao[[All,{1,2}]],f[t,a1,a2,w1,w2],{{a1,100},{a2,180},{w1,0.05},{w2,0.03}},t,Weights->1/simulacao[[All,3]]^2,VarianceEstimatorFunction->(1&)]

fitSimulacao["ParameterTable"]

chiSimul=Total[(fitSimulacao["FitResiduals"]/simulacao[[All,3]])^2 ]

Show[{ErrorListPlot[simulacao,PlotRange->All],Plot[fitSimulacao[x],{x,0,120}]}]

Um procedimento mais informativo

fixando w1 e w2, o ajuste é linear:

w1=0.05;w2=0.03

fitA=LinearModelFit[simulacao[[All,{1,2}]],{Sin[w1 t],Cos[w2 t]},t,IncludeConstantBasis->False,Weights->1/simulacao[[All,3]]^2,VarianceEstimatorFunction->(1&)]

Total[(fitA["FitResiduals"]/simulacao[[All,3]])^2 ]

soma dos quadrados dos resíduos em função do par de valores  w1 e w2 e os gráficos relacionados

q[simulacao_,w1_,w2_]:=Module[{fitA},fitA=LinearModelFit[simulacao[[All,{1,2}]],{Sin[w1 t],Cos[w2 t]},t,IncludeConstantBasis->False,Weights->1/simulacao[[All,3]]^2,VarianceEstimatorFunction->(1&)];
Total[(fitA["FitResiduals"]/simulacao[[All,3]])^2 ]
]

simulacao da medição com poucos dados

ContourPlot[q[simulacao,w1,w2],{w1,0.01,0.08},{w2,0.01,0.08}]

ContourPlot[q[simulacao,w1,w2],{w1,0.01,0.08},{w2,0.01,0.06},Contours->Table[chiSimul+n^2,{n,1,5}]]

fitSimulacao["ParameterTable"]