Confiabilidade de uma curva ajustada?

Gostaria de estimar a incerteza ou a confiabilidade de uma curva ajustada. Intencionalmente, não cito uma quantidade matemática precisa que estou procurando, pois não sei o que é.

Aqui (energia) é a variável dependente (resposta) e (volume) é a variável independente. Gostaria de encontrar a curva Energia-Volume, , de algum material. Então, eu fiz alguns cálculos com um programa de computador de química quântica para obter energia para alguns volumes de amostra (círculos verdes no gráfico). $E$ $V$ $E(V)$

Em seguida, ajustei esses exemplos de dados com a função Birch – Murnaghan : que depende de quatro parâmetros: . Também suponho que essa é a função de ajuste correta, portanto todos os erros vêm apenas do ruído das amostras. No que se segue, a função de ajuste vai ser escrito como uma função de .

E (E | V) = E_{0} + \frac{9 V_{0} B_{0}}{16} {{[{(\frac{V_{0}}{V})}^{\frac{2}{3}} - 1]}^{3} B_{0}^{'} + {[{(\frac{V_{0}}{V})}^{\frac{2}{3}} - 1]}^{2} [6 - 4 {(\frac{V_{0}}{V})}^{\frac{2}{3}}]},

$\mathbb{E}(E|V) = E_0 + \frac{9V_0B_0}{16} \left\{ \left[\left(\frac{V_0}{V}\right)^\frac{2}{3}-1\right]^3B_0^\prime + \left[\left(\frac{V_0}{V}\right)^\frac{2}{3}-1\right]^2 \left[6-4\left(\frac{V_0}{V}\right)^\frac{2}{3}\right]\right\}\;,$

E_{0}, V_{0}, B_{0}, B_{0}^{'}

$E_0, V_0, B_0, B_0'$

(\hat{E})

$(\hat{E})$

V

$V$

Aqui você pode ver o resultado (ajustando-se ao algoritmo de mínimos quadrados). A variável eixo y é e a variável do eixo-x é . A linha azul é o ajuste e os círculos verdes são os pontos de amostra. $E$ $V$

Agora preciso de alguma medida da confiabilidade (na melhor das hipóteses, dependendo do volume) dessa curva ajustada, , porque preciso calcular outras quantidades, como pressões de transição ou entalpias. $\hat{E}(V)$

Minha intutição me diz que a curva ajustada é mais confiável no meio, então acho que a incerteza (por exemplo, faixa de incerteza) deve aumentar perto do final dos dados da amostra, como neste esboço:

No entanto, qual é o tipo de medida que estou procurando e como posso calculá-lo?

Para ser preciso, na verdade, existe apenas uma fonte de erro aqui: As amostras calculadas são barulhentas devido a limites computacionais. Portanto, se eu calcular um conjunto denso de amostras de dados, elas formarão uma curva irregular.

Minha idéia para encontrar a estimativa de incerteza desejada é calcular o seguinte '' erro '' com base nos parâmetros conforme você o aprende na escola ( propagação da incerteza ):

Δ E (V) = \sqrt{{(\frac{\partial E (V)}{\partial E_{0}} Δ E_{0})}^{2} + {(\frac{\partial E (V)}{\partial V_{0}} Δ V_{0})}^{2} + {(\frac{\partial E (V)}{\partial B_{0}} Δ B_{0})}^{2} + {(\frac{\partial E (V)}{\partial B_{0}^{'}} Δ B_{0}^{'})}^{2}}

$\Delta E(V) = \sqrt{ \left(\frac{\partial E(V)}{\partial E_0} \Delta E_0\right)^2 + \left(\frac{\partial E(V)}{\partial V_0} \Delta V_0\right)^2 + \left(\frac{\partial E(V)}{\partial B_0} \Delta B_0\right)^2 + \left(\frac{\partial E(V)}{\partial B_0'} \Delta B_0'\right)^2}$ O e são fornecidos pelo software de ajuste.

Δ E_{0}, Δ V_{0}, Δ B_{0}

$\Delta E_0, \Delta V_0, \Delta B_0$

Δ B_{0}^{'}

$\Delta B_0'$

Essa é uma abordagem aceitável ou estou fazendo errado?

PS: Eu sei que também poderia resumir os quadrados dos resíduos entre minhas amostras de dados e a curva para obter algum tipo de "erro padrão", mas isso não depende do volume.

error fitting reliability uncertainty Tomilho
fonte

nenhum dos seus parâmetros é um expoente, o que é bom. Qual software NLS você usou? A maioria retornará uma estimativa para a incerteza paramétrica (que pode ser completamente irrealista se seus parâmetros forem expoentes, mas esse não for o seu caso).

DeltaIV 24/11

Não há A no lado direito da sua equação, mas ele aparece no seu gráfico. Quando você diz "quatro parâmetros", você quer dizer parâmetros no sentido estatístico (nesse caso, onde estão seus IVs) ou quer dizer variáveis (nesse caso, onde estão seus parâmetros)? Por favor, esclareça os papéis dos símbolos - o que é medido e o que são desconhecidos?

Glen_b -Reinstala Monica 25/11

Eu acho que o V é A ^ 3. foi isso que eu usei e meu enredo parecia idêntico ao dele.

Dave fournier #:

@Glen_b Acabei de assumir que o eixo Y é E na função Birch – Murnaghan enquanto o eixo x é V. Os quatro parâmetros são os quatro parâmetros na função Birch – Murnaghan. Se você assumir que recebe algo que se parece com o que ele tem.

Dave fournier:

Ah, espere, finalmente entendi. não é um operador de expectativa (como eu esperaria ver no LHS de uma equação sem um termo de erro no RHS), é a variável de resposta escrita como uma função na forma . DICA GRANDE para todos: não mostre uma equação com à esquerda de uma equação de regressão para um estatístico sem definir cuidadosamente o que você quer dizer, porque eles provavelmente assumirão que é uma expectativa.

E ()

$E()$

E

$E$

y (x)

$y(x)$

E ()

$E()$

Glen_b -Reinstate Monica

Respostas:

Este é um problema comum de mínimos quadrados!

Definindo

x = V^{- 2 / 3}, w = V_{0}^{1 / 3},

$x = V^{-2/3}, \ w = V_0^{1/3},$

o modelo pode ser reescrito

E (E | V) = β_{0} + β_{1} x + β_{2} x^{2} + β_{3} x^{3}

$\mathbb{E}(E|V) = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3$

onde os coeficientes estão algebricamente relacionados aos coeficientes originais via $\beta=(\beta_i)^\prime$

16 β = (\begin{matrix} 16 E_{0} + 54 B_{0} w^{3} - 9 B_{0} B_{0}^{'} w^{3} \\ - 144 B_{0} w^{5} + 27 B_{0} B_{0}^{'} w^{5} \\ 126 B_{0} w^{7} - 27 B_{0} B_{0}^{'} w^{7} \\ - 36 B_{0} w^{9} + 9 B_{0} B_{0}^{'} w^{9} \end{matrix}) .

$16 \beta = \pmatrix{16 E_0 + 54 B_0 w^3 - 9 B_0 B_0^\prime w^3\\ -144 B_0 w^5 + 27 B_0 B_0^\prime w^5\\ 126 B_0 w^7 - 27 B_0 B_0^\prime w^7\\ -36 B_0 w^9 + 9 B_0 B_0^\prime w^9}.$

Este é simples de resolver algebricamente ou numericamente: escolher a solução para que , e são positivos. A única razão para fazer isso é obter estimativas de , e e para verificar se eles estão fisicamente significativa. Todas as análises do ajuste podem ser realizadas em termos de . $B_0, B_0^\prime$ $w$ $B_0, B_0^\prime, w$ $E_0$ $\beta$

Essa abordagem não é apenas muito mais simples do que o ajuste não-linear, mas também é mais precisa: a matriz de variância-covariância para retornada por um ajuste não-linear é apenas uma aproximação quadrática local à distribuição da amostra desses parâmetros, enquanto que (para erros normalmente distribuídos na medição , de qualquer maneira), os resultados do OLS não são aproximações. $(E_0, B_0, B_0^\prime, V_0)$ $E$

Intervalos de confiança, intervalos de previsão, etc. podem ser obtidos da maneira usual sem a necessidade de encontrar esses valores: calcule-os em termos das estimativas e sua matriz de variância-covariância. (Até o Excel pode fazer isso!) Aqui está um exemplo, seguido pelo código (simples) que o produziu. $\hat\beta$ R

#
# The data.
#
X <- data.frame(V=c(41, 43, 46, 48, 51, 53, 55.5, 58, 60, 62.5),
                E=c(-48.05, -48.5, -48.8, -49.03, -49.2, -49.3, -49.35, 
                    -49.34, -49.31, -49.27))
#
# OLS regression.
#
fit <- lm(E ~ I(V^(-2/3)) + I(V^(-4/3)) + I(V^(-6/3)), data=X)
summary(fit)
beta <- coef(fit)
#
# Prediction, including standard errors of prediction.
#
V0 <- seq(40, 65)
y <- predict(fit, se.fit=TRUE, newdata=data.frame(V=V0))
#
# Plot the data, the fit, and a three-SEP band.
#
plot(X$V, X$E, xlab="Volume", ylab="Energy", bty="n", xlim=c(40, 60))
polygon(c(V0, rev(V0)), c(y$fit + 3*y$se.fit, rev(y$fit - 3*y$se.fit)),
        border=NA, col="#f0f0f0")
curve(outer(x^(-2/3), 0:3, `^`) %*% beta, add=TRUE, col="Red", lwd=2)
points(X$V, X$E)

Se você estiver interessado na distribuição conjunta das estimativas de parâmetros originais, é fácil simular a partir da solução OLS: simplesmente gere realizações normais multivariadas de e converta-as nos parâmetros. Aqui está uma matriz de dispersão de 2000 dessas realizações. A forte curvilinearidade mostra por que o método Delta provavelmente apresenta resultados ruins. $\beta$

whuber
fonte

Embora seja verdade que algoritmos para ajuste de modelos lineares são muito mais estáveis numericamente do que para modelos não lineares, não é verdade que exista uma diferença na precisão dos diagnósticos, desde que o algoritmo de ajuste não linear converja. Eu verifiquei e temos a mesma soma residual de quadrados para pelo menos 4 sig figs. Além disso, a parametrização linear que você escolheu é altamente confusa, de modo que nenhum dos parâmetros é significativo de acordo com o teste t. Todos os meus são. Não é realmente um grande negócio, mas divertido e pode confundir o jovem jogador.

Dave fournier #:

Além disso, eu acho que você não respondeu à pergunta do OP desde que ela afirmou que ela queria algo como limites de confiança para o volume-função entalpia

dave Fournier

@ Dave Esses são pontos importantes, obrigado. Em um problema físico, normalmente não se preocupa com significado: a teoria implica todas as variáveis. Em vez disso, trata-se de estimar os valores. Embora ambas as abordagens devam alcançar a mesma perda mínima (soma dos quadrados dos resíduos), o OLS produz distribuições corretas para a variação amostral de seus parâmetros. A abordagem não linear não. É preciso aplicar a transformação da distribuição de para , mas o uso de covariâncias dos é apenas uma aproximação.

β

$\beta$

(E_{0}, \dots)

$(E_0,\ldots)$

({\hat{E}}_{0} \dots)

$(\hat E_0\ldots)$

whuber

Seu modelo e o meu são idênticos, independentemente da parametrização. (Estou falando do modelo OLS.) É verdade que se um parâmetro específico entrar no modelo linearmente, os desvios padrão produzirão melhores limites de confiança para esse parâmetro. os desvios padrão obtidos pelo método delta serão os mesmos, sejam eles usados para parametrizar o modelo ou resolvidos como uma variável dependente. Nesse caso, a variável dependente de interesse é a função entalpia-volume e seu método delta std dev será o mesmo, independentemente de alguém usar sua parametrização ou a minha.

Dave fournier #:

@ Dave Concordo: os modelos são os mesmos, mas com parâmetros diferentes. As vantagens da formulação OLS estendem-se ao fato de que os erros normais de iid na resposta se traduzem em uma solução exata para a distribuição dos parâmetros estimados , que é facilmente transformada em uma estimativa da distribuição quadridimensional completa das estimativas dos parâmetros originais. Embora isso possa ser feito usando a parametrização original, isso exigiria mais trabalho numérico. Além disso, as vantagens conceituais de ver o modelo original são apenas OLS disfarçadas são importantes.

\hat{β}

$\hat\beta$

whuber

Existe uma abordagem padrão para isso chamada método delta. Você forma o inverso do hessiano da probabilidade logarítmica por seus quatro parâmetros. Há um parâmetro extra para a variação dos resíduos, mas ele não desempenha um papel nesses cálculos. Então você calcula a resposta prevista para os valores desejados da variável independente e calcula seu gradiente (a derivada wrt) nesses quatro parâmetros. Chame o inverso do Hessiano e o vetor gradiente . Você forma o produto da matriz vetorial $I$ $g$

- g^{t} I g

$-g^tIg$ Isso fornece a variação estimada para essa variável dependente. Pegue a raiz quadrada para obter o desvio padrão estimado. então os limites de confiança são o valor previsto + - dois desvios padrão. Isso é coisa de probabilidade padrão. no caso especial de uma regressão não linear, você pode corrigir os graus de liberdade. Você tem 10 observações e 4 parâmetros para poder aumentar a estimativa da variação no modelo multiplicando por 10/6. Vários pacotes de software farão isso por você. Eu escrevi seu modelo no Modelo AD no AD Model Builder e o ajustei e calculei as variações (não modificadas). Eles serão um pouco diferentes dos seus, porque eu tive que adivinhar um pouco os valores.

                    estimate   std dev
10   pred_E      -4.8495e+01 7.5100e-03
11   pred_E      -4.8810e+01 7.9983e-03
12   pred_E      -4.9028e+01 7.5675e-03
13   pred_E      -4.9224e+01 6.4801e-03
14   pred_E      -4.9303e+01 6.8034e-03
15   pred_E      -4.9328e+01 7.1726e-03
16   pred_E      -4.9329e+01 7.0249e-03
17   pred_E      -4.9297e+01 7.1977e-03
18   pred_E      -4.9252e+01 1.1615e-02

Isso pode ser feito para qualquer variável dependente no AD Model Builder. Um declara uma variável no local apropriado no código como este

   sdreport_number dep

e escreve o código para avaliar a variável dependente como esta

dep=sqrt(V0-cube(Bp0)/(1+2*max(V)));

Observe que isso é avaliado para um valor da variável independente 2 vezes a maior observada no ajuste do modelo. Ajuste o modelo e obtenha o desvio padrão para essa variável dependente

19   dep          7.2535e+00 1.0980e-01

Modifiquei o programa para incluir o código para calcular os limites de confiança para a função de volume de entalpia. O arquivo de código (TPL) parece

DATA_SECTION
 init_int nobs
 init_matrix data(1,nobs,1,2)
 vector E
 vector V
 number Vmean
LOC_CALCS
 E=column(data,2);
 V=column(data,1);
 Vmean=mean(V);

PARAMETER_SECTION
 init_number E0
 init_number log_V0_coff(2)
 init_number log_B0(3)
 init_number log_Bp0(3)
 init_bounded_number a(.9,1.1)
 sdreport_number V0
 sdreport_number B0
 sdreport_number Bp0
 sdreport_vector pred_E(1,nobs)
 sdreport_vector P(1,nobs)
 sdreport_vector H(1,nobs)
 sdreport_number dep
 objective_function_value f
PROCEDURE_SECTION
  V0=exp(log_V0_coff)*Vmean;
  B0=exp(log_B0);
  Bp0=exp(log_Bp0);
  if (current_phase()<4)
  f+=square(log_V0_coff) +square(log_B0);

  dvar_vector sv=pow(V0/V,0.66666667);
  pred_E=E0 + 9*V0*B0*(cube(sv-1.0)*Bp0
    + elem_prod(square(sv-1.0),(6-4*sv)));

  dvar_vector r2=square(E-pred_E);
  dvariable vhat=sum(r2)/nobs;
  dvariable v=a*vhat;
  f=0.5*nobs*log(v)+sum(r2)/(2.0*v);

  // code to calculate the  enthalpy-volume function
  double delta=1.e-4;
  dvar_vector svp=pow(V0/(V+delta),0.66666667);
  dvar_vector svm=pow(V0/(V-delta),0.66666667);
  P = -((9*V0*B0*(cube(svp-1.0)*Bp0
      + elem_prod(square(svp-1.0),(6-4*svp))))
      -(9*V0*B0*(cube(svm-1.0)*Bp0
      + elem_prod(square(svm-1.0),(6-4*svm)))))/(2.0*delta);
  H=E+elem_prod(P,V);

dep=sqrt(V0-cube(Bp0)/(1+2*max(V)));

Em seguida, reajustei o modelo para obter os desenvolvedores padrão para as estimativas de H.

29   H           -3.9550e+01 5.9163e-01
30   H           -4.1554e+01 2.8707e-01
31   H           -4.3844e+01 1.2333e-01
32   H           -4.5212e+01 1.5011e-01
33   H           -4.6859e+01 1.5434e-01
34   H           -4.7813e+01 1.2679e-01
35   H           -4.8808e+01 1.1036e-01
36   H           -4.9626e+01 1.8374e-01
37   H           -5.0186e+01 2.8421e-01
38   H           -5.0806e+01 4.3179e-01

Eles são calculados para os valores de V observados, mas podem ser facilmente calculados para qualquer valor de V.

Foi apontado que este é realmente um modelo linear para o qual existe um código R simples para realizar a estimativa de parâmetros via OLS. Isso é muito atraente, especialmente para usuários ingênuos. No entanto, desde o trabalho de Huber há mais de trinta anos, sabemos ou devemos saber que provavelmente deve-se quase sempre substituir o OLS por uma alternativa moderadamente robusta. A razão pela qual isso não é feito rotineiramente, acredito, é que métodos robustos são inerentemente não-lineares. Deste ponto de vista, os métodos OLS simples e atraentes em R são mais uma armadilha do que um recurso. Uma vantagem da abordagem do AD Model Builder é o suporte integrado à modelagem não linear. Para alterar o código dos mínimos quadrados para uma mistura normal robusta, apenas uma linha do código precisa ser alterada. A linha

    f=0.5*nobs*log(v)+sum(r2)/(2.0*v);

é alterado para

f=0.5*nobs*log(v)
  -sum(log(0.95*exp(-0.5*r2/v) + 0.05/3.0*exp(-0.5*r2/(9.0*v))));

A quantidade de sobredispersão nos modelos é medida pelo parâmetro a. Se igual a 1,0, a variação é a mesma do modelo normal. Se houver inflação da variação por valores extremos, esperamos que a seja menor que 1,0. Para esses dados, a estimativa de a é de cerca de 0,23, de modo que a variação é de cerca de 1/4 da variação do modelo normal. A interpretação é que os valores discrepantes aumentaram a estimativa de variância em um fator de cerca de 4. O efeito disso é aumentar o tamanho dos limites de confiança dos parâmetros para o modelo OLS. Isso representa uma perda de eficiência. Para o modelo de mistura normal, os desvios padrão estimados para a função de volume de entalpia são

 29   H           -3.9777e+01 3.3845e-01
 30   H           -4.1566e+01 1.6179e-01
 31   H           -4.3688e+01 7.6799e-02
 32   H           -4.5018e+01 9.4855e-02
 33   H           -4.6684e+01 9.5829e-02
 34   H           -4.7688e+01 7.7409e-02
 35   H           -4.8772e+01 6.2781e-02
 36   H           -4.9702e+01 1.0411e-01
 37   H           -5.0362e+01 1.6380e-01
 38   H           -5.1114e+01 2.5164e-01

Vê-se que há pequenas mudanças nas estimativas pontuais, enquanto os limites de confiança foram reduzidos para cerca de 60% dos produzidos pelo OLS.

O ponto principal que quero destacar é que todos os cálculos modificados ocorrem automaticamente quando um altera a linha de código no arquivo TPL.

Dave Fournier
fonte

Para o benefício do @ tomilho, eu gostaria de observar que o "método delta" é essencialmente o mesmo procedimento que a "propagação da incerteza" com a qual ele está familiarizado e que geralmente é ensinado aos cientistas - pelo menos, eles são o mesmo procedimento quando o último é feito corretamente. Uma ressalva é que a fórmula postada na pergunta ignora as correlações entre os valores estimados dos parâmetros. Ignorar as correlações é equivalente a considerar apenas os elementos diagonais de no método delta.

I

$I$

jwimberley

Também para @thyme, a propagação de incertezas / o método delta produz apenas a incerteza em . Além disso, existe qualquer tendência / variação devido ao ruído. Suponho que você esteja fazendo previsões sobre amostras físicas, cujas quantidades de energia / entalpia / outras quantidades termodinâmicas não tenham o mesmo ruído que no seu software de simulação, mas, em ambos os casos, isso adiciona uma fonte extra de variação além da a variação em ou devido à incerteza do ajuste.

E (E ∣ V)

$\mathbb{E}(E\mid V)$

E (E ∣ V)

$\mathbb{E}(E\mid V)$

E (H ∣ V)

$\mathbb{E}(H \mid V)$

jwimberley

@jwimberley, você está basicamente dizendo que dave fourier deu a fórmula para o intervalo de confiança da média (condicional), enquanto o tomilho pode estar interessado no intervalo de previsão para uma nova observação. É fácil calcular o último para o OLS. Como você o computa neste caso?

DeltaIV 25/11/16

@DeltaIV Ainda seria fácil neste caso - se o modelo de mínimos quadrados não linear estiver correto, os resíduos do ajuste serão distribuídos como . Portanto, a variação extra no intervalo de previsão é a variação dos resíduos de ajuste. Isso está relacionado à idéia no pós-script da resposta (que não é dependente de porque o modelo de ajuste não é heterocedástico). No entanto, mais importante, no ajuste vem de limites computacionais, enquanto no mundo vem de flutuações termodinâmicas, que provavelmente não são comparáveis.

E = f (V) + ϵ

$E = f(V) + \epsilon$

E - \hat{E}

$E -\hat E$

ϵ

$\epsilon$

V

$V$

ϵ

$\epsilon$

ϵ

$\epsilon$

jwimberley

@jwimberley Eu só mostrei os limites de confiança para os valores previstos correspondentes aos valores V observados simplesmente porque eles estavam disponíveis. Eu editei minha resposta para mostrar como obter limites de confiança para qualquer variável dependente.

Dave fournier #:

A validação cruzada é uma maneira simples de estimar a confiabilidade de sua curva: https://en.wikipedia.org/wiki/Cross-validation_(statistics)

A propagação da incerteza com diferenciais parciais é ótima: você realmente conhece e . No entanto, o programa que você está usando fornece apenas erros de ajuste (?). Eles serão otimistas demais (irrealisticamente pequenos). $\Delta E_{0},\Delta V_{0},\Delta B_{0}$ $\Delta B'$

Você pode calcular o erro de validação de uma vez, deixando um de seus pontos afastados do ajuste e usando a curva ajustada para prever o valor do ponto que foi deixado de fora. Repita isso para todos os pontos, para que cada um seja deixado uma vez. Em seguida, calcule o erro de validação da sua curva final (curva ajustada com todos os pontos) como uma média dos erros de previsão.

Isso mostrará apenas o quão sensível é o seu modelo para qualquer novo ponto de dados. Por exemplo, não lhe dirá quão impreciso é o seu modelo de energia. No entanto, isso será uma estimativa de erro muito mais realista, um mero erro de ajuste.

Além disso, você pode plotar erros de previsão como uma função do volume, se desejar.

Jman
fonte