2014-04-04 7 views
2

나는 내 자신의 데이터에서 헤 시안 행렬을 얻기 위해 노력하고, 나는이 개 결과이 - 라이브러리 (numDeriv)의 코드 독일인를 사용하여"Hessian"명령과 "numericHessian"명령 사이에서 신뢰해야하는 결과는 무엇입니까?

  • 라이브러리에서 코드 numericHessian를 사용하여
  • (maxLik)

Hessian의 결과는 numericHessian의 결과에 비해 매우 작습니다.

이 경우 어떤 결과를 신뢰할 수 있습니까?

특히, 사용하는 데이터는 350000에서 1100000까지이며 총 18 개의 데이터 값을 가진 9X2 매트릭스였습니다.

일종의 표준 편차 공식을 사용했고 "numericHessian"의 결과는 2X2 행렬을 사용하여 230에서 466까지의 범위 였지만 "Hessian"의 결과는 -3.42e-18에서 1.34e-17까지였습니다. 이전보다 훨씬 적습니다.

어떤 표준 편차에 대한 올바른 계산이라고 생각하십니까?

는 다음 코드이므로 :

data=read.table("C:/file.txt", header=T); 
data <- as.matrix(data); 
library(plyr) 
library(MASS) 
w1 = tail(data/(rowSums(data)),1) 
w2 = t(w1) 
f <- function(x){ 
w1 = tail(x/(rowSums(x)),1) 
w2 = t(w1) 
r = ((w1%*%cov(cbind(x))%*%w2)^(1/2)) 
return(r) 
} 
library(maxLik); 
numericHessian(f, t0=rbind(data[1,1], data[1,2])) 
library(numDeriv); 
hessian(f, rbind(data[1,1], data[1,2]), method="Richardson") 

File.txt를 다음이다 :

1     2 

137    201 

122    342 

142    111 

171    126 

134    123 

823    876 

634    135 

541    214 

423    142 

은 "numericHessian"의 결과이다 : 다음

  [,1]  [,2] 
[1,] 0.007105427 0.007105427 
[2,] 0.007105427 0.000000000 

, "Hessian"의 결과는 다음과 같습니다.

  [,1]   [,2] 
[1,] -3.217880e-15 -1.957243e-16 
[2,] -1.957243e-16 1.334057e-16 

미리 감사드립니다.

+0

여기에 대한 세부 정보가 충분하지 않습니다. 재현 할 수있는 예를 들어 주시겠습니까? –

+0

좋은 시작이지만 재현 할 수는 없습니다. –

+0

내 대답은 아래 편집 된 버전을 보았습니까? –

답변

4

reproducible example을주지 않았지만 어쨌든 시도합니다.

library(bbmle) 
x <- 0:10 
y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8) 
d <- data.frame(x,y) 
LL <- function(ymax=15, xhalf=6) 
    -sum(stats::dpois(y, lambda=ymax/(1+x/xhalf), log=TRUE)) 
fit <- mle2(LL) 
cc <- coef(fit) 

여기 MLE에서 마이너스 로그 우도 함수의 헤센 (제 유도체의 행렬)의 유한 차분 추정치 이러한 행렬들을 반전하는 파라미터의 분산 공분산 행렬의 추정치를 제공한다.

library(numDeriv) 
hessian(LL,cc) 
##    [,1]   [,2] 
## [1,] 1.296717e-01 -1.185789e-15 
## [2,] -1.185789e-15 4.922087e+00 

library(maxLik) 
numericHessian(LL, t0=cc) 
##   [,1]  [,2] 
## [1,] 0.1278977 0.000000 
## [2,] 0.0000000 4.916956 

그래서이 상대적으로 사소한 예를 들어, numDeriv::hessianmaxLik::numericHessian 매우 유사한 결과를 제공합니다. 그래서 당신이 우리에게 보여주지 않은 것이 있거나 문제의 수치에 특별한 것이 있어야합니다. 더 이상 진행하기 위해, 우리는이 명확하게 동의하는 재현 예를주세요 ...

dat <- matrix(c(137,201,122,342,142,111, 
       171,126,134,123,823,876, 
       634,135,541,214,423,142), 
     byrow=TRUE,ncol=2) 
f <- function(x){ 
    w1 <- tail(x/(rowSums(x)),1) 
    sqrt(w1%*%cov(cbind(x))%*%t(w1)) 
} 
p <- t(dat[1,1:2,drop=FALSE]) 
f(p) ## 45.25483 
numDeriv::hessian(f,p) 
##    [,1]   [,2] 
## [1,] -3.217880e-15 -1.957243e-16 
## [2,] -1.957243e-16 1.334057e-16 
maxLik::numericHessian(f,t0=p) 
##   [,1]  [,2] 
## [1,] 0.007105427 0.007105427 
## [2,] 0.007105427 0.000000000 

확인을해야합니다. 왜 그런지 모르겠지만,이 특별한 경우에 우리는 당신이 무슨 일을하는지 분석하고 하나가 바로 어떤 볼 수 있으므로, 사용자의 입력 행렬이 단일 열이기 때문에

  • , x/rowSums(x)는 사람의 벡터이다 마지막 요소 (w1 <- tail(...,1))는 단지 1입니다.
  • 그러면 표현식이 sqrt(cov(cbind(x)))으로 줄어 듭니다. 다시 말하자면, x이 1 열 매트릭스이므로, cov()은 단지 분산이고, sqrt(cov(.))은 단지 표준 편차 또는 벡터의 표준입니다.
  • 분산은 평균으로부터의 어떤 요소의 편차의 2 차 함수이므로 표준 편차는 평균 (0을 제외하고)과의 편차에서 다소 선형 적이므로 2 차 미분 값이 0이 될 것으로 기대합니다. ,

    maxLik::numericHessian(f,t0=p,eps=1e-3) 
    ##  [,1]   [,2] 
    ## [1,] 0 0.000000e+00 
    ## [2,] 0 -7.105427e-09 
    

    결론은 numDeriv는보다 정확한 (하지만 느린) 방식을 사용한다는 것입니다 : 그래서 numDeriv::hessian처럼 정답

우리는 또한 numericHessian에 대한 eps 증가하여 확인할 수 있습니다을주고있다 찾습니다 신중하다면 numericHessian에서 합리적인 답변을 얻을 수 있습니다.

+0

이 질문에 대한 답변을 받으려면 선택 표시를 클릭하십시오. –