24小时热门版块排行榜

返回列表

y1q2w3e

新虫 (初入文坛)

应助: 1 (幼儿园)
金币: 392.5
帖子: 19
在线: 8.4小时
虫号: 3512544
注册: 2014-11-01
专业: 计算机科学的基础理论

[求助] 广义帕累托分布拟合优度检验问题

有熟悉GPD分布检验的同志进来看看，本人门外汉，现在导师的研究中需要进行GPD分布检验，遇到了一些问题需要请教一下，希望各位多多指教。
在利用Cramér-von Mises 统计量W2 和Anderson- Darling 统计量A2进行GPD分布拟合优度检验的时候，首先需要产生一个如下图所示的临界值表，以便后续进行拟合优度检验。
关于产生这个表的过程，我的理解如下：
1.对于一系列的形状参数值(如0.1,0.2,0.3,...),依次选取其中一个形状参数值(尺度参数大小取为1)，随机产生100个GPD分布随机数；
  1.1.对这n个随机数利用不同的参数估计方法(如MOM,PWM,MLE,PL-moments)估计其形状参数和尺度参数；
  1.2.利用上一步的估计参数值(形状参数、尺度参数)计算，得到一个w2值和一个a2值；
对以上过程(1.1和1.2)重复N次，得到N个w2值和a2值，进而分别计算该形状参数下的w2和a2的不同的分位数，也就是图中的一行数值。最终得到临界值表。

以下是我用R语言写的函数代码：
quantile_generate<-function(shape_1,scale=1){
  #做一万次循环，每次产生100个随机数，计算得到一万个w2和a2，并分别保存在W2和A2中，
  #最后对W2和A2进行分位数计算
  W2<-c()
  A2<-c()
  for(j in 1:10^4){
m<-rgp(100,shape_1,scale) #产生100个随机数
m.mean<-mean(m)
m.var<-var(m)
m.length<-length(m)
#MOM估计参数
shape_est<-0.5*(1-((m.mean^2)/(m.var)))
scale_est<-0.5*m.mean*(1+((m.mean)^2/(m.var)))
#计算W2和A2
gpdj<-1-(1+(shape_est*(m))/(scale_est))^(-1/shape_est)
gpdj<-gpdj[which(is.na(gpdj)==FALSE)]
gpdj<-gpdj[gpdj>0]
m.length1<-length(gpdj)
i<-1:m.length1
c<-(2*i-1)/(2*m.length1)
w2<-1/(12*m.length1)+sum((gpdj-c)^2,na.rm = TRUE)
a2<--m.length1-sum((i*2-1)/m.length1*(log(gpdj)+log(1-gpdj[m.length1-i+1])),na.rm = TRUE)

W2<-c(W2,w2)
A2<-c(A2,a2)
  }
  w2<-quantile(W2,c(0.001,0.005,0.01,0.025,0.05,0.1,0.25,0.5))  #计算分位数
  a2<-quantile(A2,c(0.001,0.005,0.01,0.025,0.05,0.1,0.25,0.5)) #计算分位数
  return(list(w2,a2))
}
在形状参数取0.1时，函数的返回结果如下：

请问各位这个结果为什么这么大(相比前两张图中的临界值)？还有为什么A2会等于0啊？
比较着急，希望熟悉GPD的人尽量帮帮忙啊，谢谢各位！！！

W2临界值表.PNG

A2临界值表.PNG

函数运行结果.PNG

回复此楼

» 猜你喜欢

1楼 2016-04-01 22:01:18

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

相关版块跳转我要订阅楼主 y1q2w3e 的主题更新

返回列表