Tweedie复合泊松Gamma的R码

Question

我发现下面的R代码符合花呢复合泊松伽玛分布。我必须将它与我的399个索赔金额相匹配。我见过下面的R代码ptweedie.series(q, power, mu, phi)和dtweedie.series(y, power, mu, phi)。但是我不能完全理解代码，在将我的数据导入R中后，如何继续？提前谢谢。

Answer 1

首先要注意的是:从上面的评论中导入数据集会产生398个索赔，而不是399个索赔。其中一个比中位数索赔高出4个数量级。所以我怀疑是打字错误。在接下来的分析中，我排除了那个样本，留下了397个样本。

快速浏览一下维基百科中有关Tweedie Distributions的条目，就会发现这实际上是一组指数分布，其区别在于power参数(R文档中的xi)。Power=1产生泊松分布，power=2产生伽马分布，power=3产生逆高斯分布，以此类推。还为非整数幂定义了Tweedie分布。参数µ是平均值，φ是与方差相关的离散度参数。

因此，根据我的理解，基本的问题是，哪种power、mu和phi的组合产生的分布最适合您的索赔数据？

评估分布是否符合样本的一种方法是Q-Q图。这将绘制样本的分位数与测试分布的分位数。如果样本是按检验分布分布的，则Q-Q图应该是一条直线。在R代码中(使用X作为样本的向量)：

summary(X)        # NOTE: max/median > 1e4 !!!
#     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
# 1.00e+03 5.50e+03 1.20e+04 5.47e+05 2.50e+04 2.08e+08 
X <- X[X<max(X)]   # remove largest value (erroneous??)
hist(X,breaks=c(seq(1,1e5,1000),Inf),xlim=c(0,100000))

​

​

library(tweedie)
qqTweedie <- function(xi,p,mu,phi) {
  names <- c("Poisson","Gamma","Inverse Gaussian","Positive Stable")
  plot(qtweedie(p,xi,mu,phi),quantile(X,probs=p),
       main=paste0("Power = ",xi," (",names[xi],")"))
  qqline(X,prob=c(0.25,0.75), col="blue", lty=2,
         distribution=function(p) qtweedie(p,xi,mu,phi))
}
p <- seq(0.02,0.98,length=100)
par(mfrow=c(2,2))
lapply(c(1:4),qqTweedie,p=p,mu=1,phi=1)

​

​

伽玛分布和逆高斯分布都解释了您的数据，最高可达40,000。伽玛分布低估了较大索赔的频率，而逆高斯分布高估了它们的频率。那么让我们试试power=2.5。

par(mfrow=c(1,1))
xi <- 2.5
plot(qtweedie(p,xi,1,1),quantile(X,probs=p),main=paste0("Power = ",xi))
qqline(X,prob=c(0.25,0.75), col="blue", lty=2,
       distribution=function(p) qtweedie(p,xi,1,1))

​

​

因此，您的索赔数据似乎遵循power=2.5的tweedie分布。下一步是在给定power=2.5的情况下估计µ和phi。这是一个二维的非线性优化问题，所以我们使用nloptr软件包。事实证明，收敛性取决于初始参数是否相对接近最佳值，因此需要进行大量的试验和错误才能使nlopt(...)收敛。

library(nloptr)
F <- function(params){ # Note: xi, Q, and p are defined external to F
  mu  <- params[1]
  phi <- params[2]
  return(sum(Q - qtweedie(p,xi,mu,phi))^2)
}
xi <- 2.5
Q <- quantile(X,p) 
opt <- nloptr(x0=c(mu=1e4,phi=.01), eval_f=F, ub=c(5e4,.1), lb = c(1,0), 
              opts = list(algorithm="NLOPT_LN_COBYLA",maxeval=1e3,print_level=1))
opt$solution
# [1] 1.884839e+04 9.735325e-03

最后，我们确认解决方案确实很好地拟合了数据。

mu  <- opt$solution[1]
phi <- opt$solution[2]
par(mfrow=c(1,1))
hist(X,breaks=c(seq(1,1e5,1000),Inf),xlim=c(0,1e5))
x <- seq(1,1e5,1e3)
lines(x,dtweedie(x,xi,mu,phi),col="red")

​

​