导航

    精算后花园

    • 注册
    • 登录
    • 搜索
    • 版块
    • 最新
    • 话题
    • 热门
    • 用户
    • 群组
    1. 主页
    2. 谢远涛
    3. 最佳
    谢
    • 资料
    • 关注
    • 粉丝
    • 主题
    • 帖子
    • 最佳
    • 群组

    谢远涛 发布的最佳帖子

    • 再保险中高阶矩计算:对数正态分布与正态分布

      2.1 对数正态分布不完全k阶矩公式

      LN<-function(mu,sigma,L,U,k){
      Lk<-(log(L)-mu)/sigma-k*sigma
      Uk<-(log(U)-mu)/sigma-k*sigma
      exp(k*mu+k^2*sigma^2/2)*(pnorm(Uk)-pnorm(Lk))
      }
      

      verify the pdf

      LN(8.5,0.8,0,Inf,0)
      

      exam 2.1

      EX=LN(8.5,0.8,0,Inf,1)
      EX1=0.75*EX
      EX2=LN(8.5,0.8,0,25000,1)+25000*LN(8.5,0.8,25000,Inf,0)
      

      exam 2.2

      VARX=LN(8.5,0.8,0,Inf,2)-LN(8.5,0.8,0,Inf,1)^2
      VARX1=0.75^2*VARX
      EX2_SQ=LN(8.5,0.8,0,25000,2)+25000^2*LN(8.5,0.8,25000,Inf,0)
      VARX2=EX2_SQ-EX2^2
      VARX2
      

      2.2 正态分布不完全k阶矩公式

      2.2.0 正态分布不完全0阶矩公式:CDF

      pnorm(Inf)-pnorm(0)
      

      2.2.1 正态分布不完全1阶矩公式

      norm<-function(mu=0,sigma=1,L,U){
      L1<-(L-mu)/sigma
      U1<-(U-mu)/sigma
      mu*{pnorm(U1)-pnorm(L1)}-sigma*{dnorm(U1)-dnorm(L1)}
      }
      norm(,,0,1e10)
      norm(,,-Inf,Inf)
      

      2.2.2 正态分布不完全2阶矩公式

      norm2<-function(mu=0,sigma=1,L,U){
      L1<-(L-mu)/sigma
      U1<-(U-mu)/sigma
      (mu^2+sigma^2)*{pnorm(U1)-pnorm(L1)}-
      sigma*{(2*mu+sigma*U1)*dnorm(U1)-(2*mu+sigma*L1)*dnorm(L1)}
      }
      norm2(,,0,Inf)
      
      发布在 R
      谢
      谢远涛
      2020年5月25日 01:33
    • 按照某个变量的值拆分数据集(为循环和批量分析服务)

      /*按照变量date的值来拆分数据集
      如date=19960131生成一个数据集,date=19960228生成一个
      现在有很多个,到20080731了
      */
      /*方法一
      思路是先把distinct date的值变为宏变量,然后用data步输出。
      反复读数据,会很慢
      */
      data a;
      do date=1991 to 2011;
      output;
      end;
      run;
      data null;
      set a end=last;
      call symputx('date'||left(n),date);
      if last then call symput ('n',n);
      run;
      %macro split;
      %do i=1 %to &n;
      data b&i;
      set a;
      where date=&&date&i;
      run;
      %end;
      %mend;
      %split;

      /方法二/
      data date;
      date=input("19960131",yymmdd8.);
      output;
      date=input("19960228",yymmdd8.);
      output;
      date=input("19960309",yymmdd8.);
      output;
      run;

      data null;
      set date;
      call execute ("data data_"||strip(put(date,yymmddn8.))||";set a;if date="||date||";run;");
      run;

      /方法三/

      data a;
      format target_day 8.;inputtargetday;
      cards;
      20120730
      20120731
      20120810
      20120817
      ;run;

      %macro segmentation();
      proc sort data = a out = work.t1(keep = target_day) nodupkey;
      by target_day;
      run;
      data null;
      set work.t1 nobs=n;
      call symputx(compress("set_"||n), put(input(target_day, yymmdd8.), yymmddn8.));
      if n = 1 then call symputx('n',n);
      run;
      data %do ii = 1 %to &n.;
      data_&&set_&ii..
      %end;;
      set a;
      select (target_day);
      %do ii = 1 %to &n.;
      when ("&&set_&ii..") output data_&&set_&ii..;
      %end;
      otherwise;
      end;
      run;
      proc sql;
      drop table work.t1;
      quit;
      %mend;
      %segmentation();

      发布在 SAS
      谢
      谢远涛
      2020年5月25日 01:24
    • 罚广义线性模型:LASSO vs 岭回归

      3 罚广义线性模型 51

      install.packages("lars")
      install.packages("ridge")
      

      (1) OLS回归

      library(ridge)
      data(GenCont)
      model.ols <- lm(Phenotypes ~ ., data = as.data.frame(GenCont))
      summary(model.ols)
      

      提取回归系数

      coef.ols <- coef(model.ols)
      

      查看不等于0的回归系数

      coef.ols[coef.ols!=0]  
      

      (2) ridge regression(岭回归)

      model.rid <- linearRidge(Phenotypes ~ ., seq(0, 10, length = 100), data = as.data.frame(GenCont))
      

      查看结果

      summary(model.rid)
      

      提取系数

      coef返回的是标准化系数; 直接输入model.rid返回的是原始数据系数
      coef.rid <- coef(model.rid)
      
      model.rid2 <- linearRidge(Phenotypes ~ ., data = as.data.frame(GenCont))
      summary(model.rid2)
      coef.rid <- coef(model.rid2)
      

      查看不等于0的系数

      coef.rid[coef.rid!=0]                           
      

      (3) 做 lasso regression

      模型设定

      library(lars)
      dt<-as.matrix(GenCont)
      X<-dt[,-1]
      Y<-dt[,1]
      model.lasso1 <- lars(X, Y)    
      model.lasso2 <- lars(X, Y, type='lasso') 
      model.lasso3 <- lars(X, Y, type='for')  # Can use abbreviations
      

      看变量选入顺序

      model.lasso3
      summary(model.lasso3)
      

      Cp的含义:衡量多重共线性,其取值越小越好,这里取到第3步使得Cp值最小

      画个图

      plot(model.lasso3)
      set.seed(12345)
      

      做交叉验证

      CV.lasso <- cv.lars(X, Y, K=10)   
      
      (best <- CV.lassoindex[which.min(CV.lassocv)])
      

      选择最好的效果

      (coef.lasso <- coef.lars(model.lasso3, mode='fraction', s=best))
      

      命名

      names(coef.lasso) <- colnames(dt)[-1]
      

      查看结果

      coef.lasso[coef.lasso!=0]
      
      发布在 R
      谢
      谢远涛
      2020年5月25日 00:55