적합 곡선의 신뢰성? \ hat {E}

적합 곡선의 불확실성 또는 신뢰성을 추정하고 싶습니다. 나는 그것이 무엇인지 모르기 때문에 내가 찾고있는 정확한 수학적 양을 의도적으로 언급하지 않습니다.

여기서 (에너지)는 종속 변수 (응답)이고 (볼륨)는 독립 변수입니다. 일부 재료 의 에너지-볼륨 곡선 를 찾고 싶습니다 . 그래서 나는 양자 화학 컴퓨터 프로그램을 사용하여 일부 샘플 볼륨 (플롯의 녹색 원)의 에너지를 얻기 위해 계산을했습니다. $E$

이자형

$V$

$E (V)$

E(V)

그런 다음 I은 이러한 데이터 샘플 끼워 자작 Murnaghan 기능 :

에 따라 다릅니다. 네 개의 매개 변수 : . 또한 이것이 올바른 피팅 함수라고 가정하므로 모든 오류는 샘플의 노이즈에서 발생합니다. 다음에, 적합 함수 는 의 함수로 작성됩니다 .

E (E | V) = E 0 + 9 V 0 B 0 16 ⎧ ⎩ ⎨ [(V 0 V) 2 3 - 1] 3 B' 0 + [(V 0 V) 2 3 - 1] 2 [6 - 4 (V 0 V) 2 3] ⎫ ⎭ ⎬,

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

E0,V0,비0,비0‘

$(\hat{E})$

(이자형^)

$V$

여기에서 결과를 볼 수 있습니다 (최소 제곱 알고리즘에 적합). y 축 변수는 $E$

이자형

이고 x 축 변수는 $V$

입니다. 파란색 선이 적합하고 녹색 원이 샘플 점입니다.

이제 전이 압력이나 엔탈피와 같은 추가 수량을 계산해야하기 때문에이 적합 곡선 의 신뢰성에 대한 측정 값이 $\hat{E} (V)$

이자형^(V)

필요합니다.

내 직감에 따르면 적합 곡선이 중간에서 가장 신뢰할 수 있다고 알려주 므로이 스케치에서와 같이 불확실성 (예 : 불확실성 범위)이 샘플 데이터의 끝 부분에서 증가해야한다고 생각합니다.

그러나 내가 찾고있는이 종류의 측정은 무엇이며 어떻게 계산할 수 있습니까?

정확히 말하자면 여기에는 실제로 하나의 오류 소스 만 있습니다. 계산 된 샘플은 계산 한계로 인해 노이즈가 발생합니다. 따라서 밀도가 높은 데이터 샘플 세트를 계산하면 울퉁불퉁 한 곡선을 형성합니다.

원하는 불확실성 추정치를 찾는 나의 생각은 학교에서 배울 때 매개 변수를 기반으로 다음과 같은“오류 ”를 계산하는 것입니다 ( 불확실성의 전파 ).

Δ E (V) = (\partial 이자형 ( V ) \partial 이자형 0 Δ E 0) 2 + (\partial 이자형 ( V ) \partial V 0 Δ V 0) 2 + (\partial 이자형 ( V ) \partial 비 0 Δ B 0) 2 + (\partial 이자형 ( V ) \partial 비 ‘ 0 Δ B ‘ 0) 2 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \sqrt

및 , 피팅 소프트웨어에 의해 주어진다. $Δ E_{0}, Δ V_{0}, Δ B_{0}$

Δ이자형0,ΔV0,Δ비0

$Δ B_{0}^{'}$

Δ비0‘

그것이 수용 가능한 접근입니까, 아니면 잘못하고 있습니까?

추신 : 나는 데이터 샘플과 곡선 사이의 잔차 제곱을 합산하여 일종의“표준 오류 ”를 얻을 수는 있지만 볼륨에 의존하지는 않습니다.

답변

이것은 일반적인 최소 제곱 문제입니다!

정의

x = V - 2 / 3, w = V 1 / 3 0,

모델을 다시 쓸 수 있습니다

E (E | V) = β 0 + β 1 x + β 2 엑스 2 + β 삼 엑스 삼

계수 는 다음을 통해 원래 계수와 대수적으로 관련됩니다. $β = (β_{i})^{'}$

β=(β나는)‘

(16) β = ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ 16 E 0 + 54 B 0 승 삼 - 9 B 0 비 ‘ 0 승 삼 - 144 B 0 승 5 + 27 B 0 비 ‘ 0 승 5 126 B 0 승 7 - 27 B 0 비 ‘ 0 승 7 - 36 B 0 승 9 + 9 B 0 비 ‘ 0 승 9 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ .

이것은 대수적으로나 수치 적으로 풀기 쉽습니다 : , 가 양수인 해를 선택하십시오. 이를 수행하는 유일한 이유는 및 추정값을 획득 하고 물리적으로 의미가 있는지 확인하기위한 것입니다. 적합도에 대한 모든 분석은 로 수행 할 수 있습니다 . $B_{0}, B_{0}^{'}$

비0,비0‘

$w$

승

$B_{0}, B_{0}^{'}, w$

비0,비0‘,승

$E_{0}$

이자형0

$β$

위한 분산 공분산 행렬이 접근법은 또한보다 정확한 만 훨씬 간단한 비선형 피팅 이하이다 비선형 맞춤에 의해 반환이에만 로컬 차 근사 이러한 파라미터의 분포를 샘플링하는 반면 ( 측정에서 정상적으로 분포 된 오차의 경우) OLS 결과는 근사치가 아닙니다. $(E_{0}, B_{0}, B_{0}^{'}, V_{0})$

(이자형0,비0,비0‘,V0)

$E$

이자형

$\hat{β}$

β^

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)

$β$

답변

$I$

나는

$g$

지

- g 티 나는 지

그러면 해당 종속 변수에 대한 추정 분산이 제공됩니다. 추정 된 표준 편차를 얻으려면 제곱근을 취하십시오. 신뢰 한계는 예측 된 값 +-2 표준 편차입니다. 이것은 표준 가능성입니다. 비선형 회귀의 특수한 경우 자유도를 교정 할 수 있습니다. 10 개의 관측치와 4 개의 매개 변수가 있으므로 10/6을 곱하여 모형의 분산 추정치를 증가시킬 수 있습니다. 여러 소프트웨어 패키지가이를 수행합니다. AD Model Builder의 AD Model로 모델을 작성하고 적합하게하고 수정되지 않은 분산을 계산했습니다. 값을 조금 추측해야하기 때문에 그들은 당신과 약간 다를 것입니다.

                    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

이는 AD Model Builder의 모든 종속 변수에 대해 수행 할 수 있습니다. 다음과 같이 코드에서 적절한 지점에 변수를 선언합니다.

   sdreport_number dep

다음과 같이 종속 변수를 평가하는 코드를 작성합니다.

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

이는 모형 피팅에서 관찰 된 최대 변수의 2 배인 독립 변수 값에 대해 평가됩니다. 모형을 적합하면이 종속 변수에 대한 표준 편차를 얻습니다.

19   dep          7.2535e+00 1.0980e-01

엔탈피-볼륨 함수의 신뢰 한계를 계산하기위한 코드를 포함하도록 프로그램을 수정했습니다. 코드 (TPL) 파일은 다음과 같습니다.

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)));

그런 다음 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

이 값은 관측 된 V 값에 대해 계산되지만 모든 V 값에 대해 쉽게 계산할 수 있습니다.

이것은 실제로 OLS를 통해 파라미터 추정을 수행하는 간단한 R 코드가있는 선형 모델이라는 것이 지적되었습니다. 이것은 특히 순진한 사용자에게 매우 매력적입니다. 그러나 30 년 전에 Huber의 작업을 통해 OLS를 거의 항상 대체 할만한 대안으로 대체해야한다는 것을 알고 있거나 알아야합니다. 이것이 일상적으로 이루어지지 않는 이유는 강력한 방법이 본질적으로 비선형이기 때문입니다. 이 관점에서 R의 간단한 매력적인 OLS 방법은 기능이 아니라 함정에 가깝습니다. AD Model Builder 접근법의 advabtage는 비선형 모델링에 대한 지원을 기본적으로 제공합니다. 최소 제곱 코드를 강력한 일반 혼합으로 변경하려면 한 줄의 코드 만 변경하면됩니다. 라인

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

로 변경되었습니다

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에 의해 측정된다. 1.0과 같으면 분산은 일반 모형과 동일합니다. 특이 치에 의한 분산의 인플레이션이있는 경우 a가 1.0보다 작을 것으로 예상합니다. 이 데이터에서 a의 추정치는 약 0.23이므로 분산은 정규 모형의 분산의 약 1/4입니다. 해석은 특이 치가 분산 추정치를 약 4 배 증가 시켰다는 것입니다. 이것의 결과는 OLS 모델의 모수에 대한 신뢰 한계의 크기를 증가시키는 것입니다. 이것은 효율성의 손실을 나타냅니다. 정규 혼합 모형의 경우 엔탈피-볼륨 함수의 추정 표준 편차는 다음과 같습니다.

 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

포인트 추정치에는 약간의 변화가 있지만, 신뢰 한계는 OLS에서 생성 된 것의 약 60 %로 줄었습니다.

내가 만들고 싶은 요점은 TPL 파일에서 한 줄의 코드를 변경하면 수정 된 모든 계산이 자동으로 수행된다는 것입니다.

답변

교차 유효성 검사는 곡선의 안정성 을 평가 하는 간단한 방법입니다 .
https://en.wikipedia.org/wiki/Cross-validation_(statistics)

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

Δ이자형0,ΔV0,Δ비0

$Δ B^{'}$

Δ비‘

포인트 중 하나를 피팅에서 멀리두고 피팅 곡선을 사용하여 누락 된 포인트의 값을 예측하여 1 배 유효성 검사 오류를 계산할 수 있습니다. 모든 포인트에 대해이 작업을 반복하여 각 포인트가 한 번만 사라지도록하십시오. 그런 다음 최종 곡선의 유효성 검사 오류 (모든 점에 맞는 커브)를 평균 예측 오류로 계산하십시오.

이것은 새로운 데이터 포인트에 대한 모델의 민감도를 알려줍니다. 예를 들어, 에너지 모델이 얼마나 부 정확한지 알려주지 않습니다. 그러나 이것은 단순한 피팅 오차보다 훨씬 현실적인 오차 추정치입니다.

또한 원하는 경우 예측 오류를 부피의 함수로 플로팅 할 수 있습니다.

답변

Exit mobile version

How IT

언제든지 물어보세요.

적합 곡선의 신뢰성? \ hat {E}

답변

답변

답변

답변