0

我想拟合变换帕累托分布，然后需要计算以下数据的 Hessian 矩阵。

library(stats4)
library(MASS)
library(vcd)         # for goodness of fit test
library(pracma)      # for hessain matrix
library(numDeriv)   

# Data from Exceedances of Wheaton River flood data.

x =c(1.7, 2.2,   14.4, 1.1,  .4, 20.6, 5.3, .7, 1.9, 13, 12, 9.3, 1.4, 18.7, 8.5, 25.5, 11.6, 
     14.1, 22.1, 1.1, 2.5, 14.4, 1.7, 37.6 ,.6, 2.2, 39, .3, 15, 11, 7.3, 22.9, 1.7, .1, 1.1, 
     .6, 9, 1.7, 7,  20.1, .4, 2.8, 14.1, 9.9, 10.4, 10.7, 30, 3.6, 5.6, 30.8, 13.3, 4.2, 25.5, 
     3.4, 11.9, 21.5, 27.6, 36.4, 2.7, 64, 1.5, 2.5, 27.4, 1, 27.1, 20.2, 16.8, 5.3, 9.7, 27.5, 
     2.5, 27)



k=.35                       # guessed vales
gamma=.1                     # minimum vales of x ,p, 0
lambda=-.95                   # guessed vale
theta=c(k,lambda)           
fn=function(k,lambda) 
{
  n=length(x)
  -n*log(k)-n*(k)*log(.1)+(k+1)*sum(log(x))-sum(log((1-lambda)+2*lambda*((.1/x)^(x))))
}
result=nlm(fn, p=c(1), theta,  hessian=TRUE, print.level=2 ) # minimization
print(result)
result1=solve(result$hessian)  # inverse of Hesssain approx
print(result1)

此代码的结果仅提供一个不正确的值，我也需要 2 x 2 hessian 矩阵。提前致谢。

4

1 回答 1

1

fn应该有一个长度为 2 的参数，并且nlmargs 需要固定。由于我们正在记录，k我们添加了一个条件来防止k下降到接近零或更低。

fn=function(p) 
{
  k <- p[1]
  if (k < 1e-10) return(10^10) # optional: will eliminate the warnings
  lambda <- p[2]
  n=length(x)
  -n*log(k)-n*(k)*log(.1)+(k+1)*sum(log(x))-sum(log((1-lambda)+2*lambda*((.1/x)^(x))))
}

result=nlm(fn, theta,  hessian=TRUE, print.level=2 ) # minimization

result1=solve(result$hessian)  # inverse of Hesssain approx
print(result1)

给予：

             [,1]         [,2]
[1,] 0.0008266354 0.0000000000
[2,] 0.0000000000 0.0009375602
于 2016-11-21T00:36:16.003 回答