拟合转换后的Pareto分布和Hessian矩阵的计算

我想拟合转换后的Pareto分布，然后需要计算以下数据的Hessian矩阵。拟合转换后的Pareto分布和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乘2 Hessian矩阵。在此先感谢。

来源

2016-11-20 Iftikhar

fn应该有一个长度为2的参数，并且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

来源

2016-11-21 00:36:16

拟合转换后的Pareto分布和Hessian矩阵的计算

回答

相关问题