Quão significativos ou úteis são os erros de parâmetro produzidos ao executar um LinearModelFit ou NonlinearModelFit não ponderado?
Esta pode ser uma questão que oscila à beira de pertencer mais ao reino das estatísticas e SE com validação cruzada , mas também estou especificamente interessado nas rotinas de adaptação do Mathematica.
Normalmente, se eu quiser ajustar um modelo a alguns dados usando um NonlinearModelFit
ou outro LinearModelFit
, terei alguns erros associados ao meu$y$-dados que uso para pesar os ajustes. Esses erros podem ser simplesmente o erro padrão adquirido de medições repetidas, ou posso saber algo sobre os processos físicos e posso atribuir pesos.
Por exemplo, Weights->1/YDataErrors^2
e sempre defino meu estimador de variância como VarianceEstimatorFunction -> (1 &)
. Posso obter meus erros de parâmetro da matriz de covariância ou simplesmente com MyFit["ParameterErrors"]
.
No entanto, em alguns casos, pode não haver erros para os dados que deseja ajustar, o que significa que não se pode fornecer pesos da maneira que descrevi acima. Minha pergunta é, então, quão confiáveis - ou mais importante - quão fisicamente / estatisticamente significativos são os erros de parâmetro para um ajuste não ponderado no Mathematica?
Respostas
Por exemplo, se houver duas fontes de erro, digamos um erro de medição e um erro de falta de ajuste, usar os pesos com base nos erros de medição pode resultar em subestimações grosseiras dos erros padrão. Considere o seguinte modelo:
$$y=a+b x +\gamma + \epsilon$$
Onde $y$ é a resposta medida, $x$ é o preditor, $a$ e $b$ são constantes a serem estimadas, $\gamma$ é o erro de medição repetido com $\gamma \sim N(0,\sigma_{ME})$e $\epsilon$ é o erro de falta de ajuste com $\epsilon \sim N(0,\sigma)$ e todos os erros são considerados independentes.
Primeiro defina alguns parâmetros específicos:
(* Measurement error standard deviation *)
σME = 10;
(* Lack-of-fit error standard deviation *)
σ = 20;
(* Regression coefficients *)
a = 1;
b = 1;
Gere e plote alguns dados:
n = 100;
x = Range[n];
SeedRandom[12345];
measurementError = RandomVariate[NormalDistribution[0, σME], n];
lackOfFitError = RandomVariate[NormalDistribution[0, σ], n];
y = a + b x + measurementError + lackOfFitError;
data = Transpose[{x, y}];
data2 = {#[[1]], Around[#[[2]], σME]} & /@ data;
ListPlot[data2]
Agora, considere dois modelos lineares diferentes onde lm1
está o que você sugere e lm2
é o que eu sugiro:
lm1 = LinearModelFit[data, z, z, Weights -> 1/ConstantArray[σME^2, n],
VarianceEstimatorFunction -> (1 &)];
lm2 = LinearModelFit[data, z, z];
lm1["ParameterTable"]
lm2["ParameterTable"]
As estimativas dos parâmetros são idênticas, mas os erros padrão para lm1
são menos da metade do tamanho daqueles para lm2
. Qual está correto?
A matriz de covariância "verdadeira" dos estimadores de mínimos quadrados de a
e b
para este modelo é
$$\left(\sigma ^2+\sigma_{ME}^2\right) \left(X^T.X\right)^{-1}$$
Onde $X$é a matriz de design. No código do Mathematica , o erro padrão para b
é
X = Transpose[{ConstantArray[1, n], Range[n]}]
Sqrt[(σME^2 + σ^2) Inverse[Transpose[X].X][[2, 2]]] // N
(* 0.0774635 *)
Isso combina muito bem com lm2
.
Este é um exemplo ligeiramente artificial em que todos os erros padrão de medição são idênticos porque as funções de regressão do Mathematica permitem apenas um único termo de erro. E tendo os erros padrão de medição idênticos, isso resulta em um modelo equivalente com um único erro.
No entanto, mesmo quando os desvios-padrão de medição variam consideravelmente, a questão sobre a ponderação inadequada de forma que não corresponda à estrutura de erro do modelo permanece.
As rotinas de regressão do Mathematica ainda não são adequadas para modelos com mais de uma fonte de erro. Eu gostaria que fossem.