継承されたパレート分布とヘッセ行列の計算

Question

私はTransmuted Pareto Distributionに適合したいので、計算する必要があります 次のデータのヘッセ行列。

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) 

このコードの結果は、正確ではない1つの値を提供します.2ヘッセ行列は2つ必要です。 ありがとうございます。

Answer 1

fnは、長さ2の引数を1つ持ち、nlm引数を固定する必要があります。 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