| 查看: 4060 | 回复: 0 | |||
[求助]
广义帕累托分布拟合优度检验问题
|
|
有熟悉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 |
回复此楼» 猜你喜欢
filecode=后面第一个是大写字母
已经有10人回复
-
售SCI一区T0P文章，我:8O.55.1.O.54,科目全,可伽急
已经有5人回复
-
售SCI一区文章，我:8.O.551.O.5.4,科目全,可伽急
已经有3人回复
-
欢迎发来filecode的Mz6后的代码验证其规律
已经有81人回复
-
售SCI一区文章，我:8.O.551.O.5.4,科目全,可伽急
已经有7人回复
-
感觉是下周放榜了
已经有7人回复
-
时间戳变了，能看出什么问题？
已经有4人回复
-
哪位老哥知道今年的国自然具体哪一天放榜？
已经有15人回复
-
售SCI文章，我:8O5.5.1.O.54,科目齐全,可+急
已经有6人回复
-
售一区SCI文章T0P，我:8O.551.O54,科目全,可十急
已经有4人回复
25