用于Tweedie复合Poisson Gamma的R代码

Question

我发现以下R代码适合Tweedie复合Poisson Gamma分布.我必须适应我的399索赔金额.我见过以下R代码ptweedie.series(q, power, mu, phi)和dtweedie.series(y, power, mu, phi).但是,我无法完全理解代码,并在将数据导入R后,如何继续？提前致谢.

Answer 1

首先注意:从上面的评论中导入数据集产生了398件索赔,而不是399件.其中一件比中位数索赔大4个数量级.所以我怀疑是一个错字.在随后的分析中,我排除了样本,留下397.

快速浏览一下Tweedie Distributions的维基百科条目,可以看出这实际上是一个由power参数区分的指数分布族(xi在R文档中).功率= 1产生泊松分布,功率= 2产生伽马分布,功率= 3产生逆高斯分布,等等.Tweedie分布也定义为非整数幂.参数mu是平均值,phi是与方差相关的色散参数.

因此,根据我的理解,基本问题是权力,亩和phi的哪种组合产生最适合您的索赔数据的分布？

评估分布是否适合样本的一种方法是QQ图.这绘制了样本的分位数与测试分布的分位数.如果样本作为测试分布分布,则QQ图应为直线.在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)

Gamma和反向高斯分布都可以解释您的数据,最多可达到约40,000.Gamma分布低估了较大声明的频率,而逆高斯分布则高估了它们的频率.那么让我们试试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))

因此,您的声明数据似乎遵循功率= 2.5的tweedie发行版.下一步是估算mu和phi,给定功率= 2.5.这是2维的非线性优化问题,因此我们使用包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")