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

Vito, 7/6/2016

Needs["ErrorBarPlots`"]

dados experimentais

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

graficoDados=ErrorListPlot[seno,PlotRange->All]

Ajuste dos parâmetros

fitSemCuidado=NonlinearModelFit[seno[[All,{1,2}]],{α Sin[ω t]},{α,ω},{t}]

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

fitSemPesos=NonlinearModelFit[seno[[All,{1,2}]],{α Sin[ω t]},{{α,1000},{ω,0.06}},{t}]

fitSemPesos["ParameterTable"]

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

fitCerto=NonlinearModelFit[seno[[All,{1,2}]],{α Sin[ω t]},{{α,1000},{ω,0.06}},{t},Weights->1/seno[[All,3]]^2,VarianceEstimatorFunction->(1&)]

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

simulacaoMedicao[seno,1000,0.06]

fit=NonlinearModelFit[simulacaoMedicao[seno,1267,0.0615],{α Sin[ω t]},{{α,1000},{ω,0.06}},{t},Weights->1/seno[[All,3]]^2,VarianceEstimatorFunction->(1&)]

fit["ParameterTable"]

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

fit["ParameterTableEntries"]

fit["ParameterTableEntries"][[All,1]]

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.

repeticoes=Table[simulaFit[seno, 1268.,0.0615],{10000}]

GraphicsRow[ {Histogram[repeticoes[[All,1]] ],Histogram[repeticoes[[All,2]] ]} ]

Mean[repeticoes]

StandardDeviation[repeticoes]/100

(Mean[repeticoes]-{1268.,0.0615})/(StandardDeviation[repeticoes]/100)

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:

StandardDeviation[repeticoes]

fitCerto["ParameterTable"]

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

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

primeiro, um exemplo no caso de um único parâmetro não-linear

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

Plot[qChapeu[seno,ω],{ω,0.05,0.07},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]}]

Não linear em dois dos parâmetros

forma funcional

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

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

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 ]

Clear[w1,w2]

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