用mle分段拟合损失分布,并计算模拟随机数



  • 因为保险损失程度分布是long tail的,所以需要有时需要对不同段的损失拟合不同分布。
    这里的思路是先分段,假设好损失出现在每段的概率,然后对每段数据进行分别拟合。
    分别拟合会遇到Density的积分加起来不等于1的情况,因此,还需要对每段的Density函数用一个系数进行调整,来达到把总和调整为1的效果。



  • 样本代码如下,假设分三段 [0-50], [50-120], [120-250]; 每段的概率是0.2, 0.5, 0.3;

    FitThreeDBs这个函数会在每个模型跑出来的参数里面加一个系数,这个系数要×在Density Function上,来实现把Density总和调整到1的效果。

    # 四个参数,数据,分段点,拟合的分布名,以及每段的概率
    FitThreeDBs = function(data,breaks,fittedDBs,probabilities){
    #把拟合数据分段
      data2cut<- 
        cut(data, breaks = breaks,
            include.lowest = TRUE)
    #储存每段的模型
      model = list()
      for(i in 1:(length(breaks)-1)){
        if(length(data)<2){
    #如果某段数据量小于2,那么不输出模型
          model= c(model,"")
        print("warning")
        }
        else{
    #拟合这段数据的模型
          temp_model = fitdistr(unlist(split(data,data2cut)[i]), fittedDBs[i])
          
        }
    #计算调整系数,调整系数应该等于在拟合数据范围我们想要的概率,除以拟合出的分布在拟合数据范围的概率
        point1 = do.call(paste0("p",fittedDBs[i]), list(breaks[i+1],temp_model$estimate[1],temp_model\$estimate[2]))
        point2 =  do.call(paste0("p",fittedDBs[i]), list(breaks[i],temp_model$estimate[1],temp_model\$estimate[2])) 
        temp_model["coefficient"] <- probabilities[i]/(point1-point2)
        model[[length(model)+1]]= temp_model
      }
      return(model)
    }
    
    #Simulation data是拟合数据
    modelCollections = FitThreeDBs(simulation_data,c(0,50,120,250),c("gamma","weibull","gamma"),c(0.2,0.5,0.3))
    
    
    


  • 验证了一下,总和确实为1

    sum = pgamma(50,shape=modelCollections[[1]]\$estimate["shape"], rate=modelCollections[[1]]\$estimate["rate"])*modelCollections[[1]]['coefficient'][[1]]
    sum = sum + (pweibull(120,modelCollections[[2]]\$estimate["shape"],modelCollections[[2]]\$estimate["scale"])- 
                          pweibull(50,modelCollections[[2]]\$estimate["shape"]
                                   ,modelCollections[[2]]\$estimate["scale"]))*modelCollections[[2]]['coefficient'][[1]]
    sum = sum + (pgamma(250,modelCollections[[3]]$estimate["shape"]
                        ,modelCollections[[3]]\$estimate["rate"])- pgamma(120,modelCollections[[3]]\$estimate["shape"]
                                                                         ,modelCollections[[3]]\$estimate["rate"]))*modelCollections[[3]]['coefficient'][[1]]
    
    #输出Sum,结果为1
    


  • 那么,如何用拟合出来的分布做Simulation呢,就要用到一个知识点,那就是可以先生成服从0到1之间均匀分布的随机数,再Map到拟合分段函数的CDF $F(x)$上,然后计算该概率下的quantile就是最后的随机数。

    分段函数的CDF $F(x)$是由每段函数原本的CDF通过改变截距和形状来的,那么我们就可以把它还原到原来函数的CDF,用R的功能计算出quantile来

    #生成1000个符合0-1之间均匀分布的随机数
    init_numebrs = runif(1000)
    
    simulation = function(init){
    #这两个break就是CDF的分段点
      break1=  0.2
      break2 = 0.2+0.5
     #这里就是把改变形状截距的CDF还原回去,然后计算quantile
      if(init<break1){
        output = qgamma(init/modelCollections[[1]]['coefficient'][[1]],shape=modelCollections[[1]]\$estimate["shape"], rate=modelCollections[[1]]\$estimate["rate"])
      }else if(init>break2){
        output = qgamma((init-break2)/modelCollections[[3]]['coefficient'][[1]]+pgamma(120,modelCollections[[3]]\$estimate["shape"],modelCollections[[3]]\$estimate["rate"]),shape=modelCollections[[3]]\$estimate["shape"], rate=modelCollections[[3]]\$estimate["rate"])
      }else{
        output = qweibull((init-break1)/modelCollections[[2]]['coefficient'][[1]]+pweibull(50,modelCollections[[2]]\$estimate["shape"],modelCollections[[2]]\$estimate["scale"]),shape=modelCollections[[2]]\$estimate["shape"], scale=modelCollections[[2]]\$estimate["scale"])
        
      }
    }
    
    simulation_result=lapply(init_numebrs,simulation)
    
    dt = data.frame("data"=unlist(simulation_result))
    p = ggplot(data = data.frame(dt), aes(x = data)) + geom_histogram(binwidth=10,fill="#69b3a2", color="#e9ecef", alpha=0.3)
    

    拟合出来的分布大概是以下形状
    fce1ef67-6829-4088-8cf6-8153b2f52b3e-image.png


登录后回复