Quão significativos ou úteis são os erros de parâmetro produzidos ao executar um LinearModelFit ou NonlinearModelFit não ponderado?

Aug 16 2020

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 NonlinearModelFitou 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^2e 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

4 JimB Aug 16 2020 at 19:20

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 lm1está 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 lm1sã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 ae bpara 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.