Tweaked simulation code.

Makefile

@@ -50,4 +50,4 @@ shiny:
 	R -e 'shiny::runApp("scripts/shiny.r")'
 
 test:
-	R -e 'source("patternBreak/tests/run_tests.r")'
+	R -e 'source("claimsBoot/tests/run_tests.r")'

claimsBoot/R/RcppExports.R

@@ -18,7 +18,7 @@ NULL
 #'
 #' @param triangle Incremental claims triangle
 #' @export
-glmBoot <- function(triangle, n_boot, boot_type, opt, seed = 42L) {
+glmBoot <- function(triangle, n_boot, boot_type, opt = "null", seed = 42L) {
     .Call('_claimsBoot_glmBoot', PACKAGE = 'claimsBoot', triangle, n_boot, boot_type, opt, seed)
 }
 

claimsBoot/R/utils.r

@@ -8,7 +8,7 @@ mackConfig <- function(ndev,
   if (boot.type != "pairs") {
     conds <- c(TRUE, FALSE)
     if (boot.type == "residuals") {
-      opt <- as.double(key[c("standardised", "studentised")])
+      opt <- as.double(key[c("standardised", "studentised", "lognormal")])
     } else if (boot.type == "parametric") {
       opt <- as.double(key[c("normal", "gamma")])
     }
@@ -115,7 +115,7 @@ mackPost <- function(res.list, res.names, sim.type) {
       "cond",
       "reserve"
     )
-  } else {
+  } else if (sim.type == "calendar") {
     out.names <- c(
       "boot.type",
       "outlier.diagidx",
@@ -126,6 +126,17 @@ mackPost <- function(res.list, res.names, sim.type) {
       "cond",
       "reserve"
     )
+  } else {
+    out.names <- c(
+      "boot.type",
+      "outlier.rowidx",
+      "mean.factor",
+      "sd.factor",
+      "excl.rowidx",
+      "opt",
+      "cond",
+      "reserve"
+    )
   }
 
   if ("pairs" %in% res.names) {
@@ -289,7 +300,7 @@ glmPost <- function(res.list, res.names, sim.type) {
       "opt",
       "reserve"
     )
-  } else {
+  } else if (sim.type == "calendar") {
     out.names <- c(
       "boot.type",
       "outlier.diagidx",
@@ -298,6 +309,15 @@ glmPost <- function(res.list, res.names, sim.type) {
       "opt",
       "reserve"
     )
+  } else {
+    out.names <- c(
+      "boot.type",
+      "outlier.rowidx",
+      "factor",
+      "excl.rowidx",
+      "opt",
+      "reserve"
+    )
   }
 
   if ("residuals" %in% res.names) {

claimsBoot/inst/unit_tests/_runit.glm.r

@@ -5,7 +5,7 @@ test.glmBoot <- function() {
 
 test.glmSim <- function() {
     ndev <- dim(test.triangle)[1]
-    config <- patternBreak:::glmConfig(ndev, factors = seq(0.5, 1.5, by = 0.25))
+    config <- claimsBoot:::glmConfig(ndev, factors = seq(0.5, 1.5, by = 0.25))
     test.triangle <- cum2incr(test.triangle)
     results <- glmSim(test.triangle, 1e3, config, "single")
     checkTrue(!any(is.na(results$reserve)))

claimsBoot/inst/unit_tests/_runit.rng.r

@@ -1,3 +1,3 @@
 test.rng <- function() {
-    checkEquals(patternBreak:::validate_rng(1e3), "Success")
+    checkEquals(claimsBoot:::validate_rng(1e3), "Success")
 }

claimsBoot/src/constants.h

@@ -14,9 +14,8 @@ inline const int RESID = 2;
 inline const int PAIRS = 3;
 
 inline const int STANDARDISED = 1;
-inline const int MODIFIED = 2;
-inline const int STUDENTISED = 3;
-inline const int LOGNORMAL = 4;
+inline const int STUDENTISED = 2;
+inline const int LOGNORMAL = 3;
 
 inline const int SINGLE = 1;
 inline const int CALENDAR = 2;
@@ -30,7 +29,6 @@ inline Rcpp::List key = Rcpp::List::create(
   Rcpp::Named("residuals") = RESID,
   Rcpp::Named("pairs") = PAIRS,
   Rcpp::Named("standardised") = STANDARDISED,
-  Rcpp::Named("modified") = MODIFIED,
   Rcpp::Named("studentised") = STUDENTISED,
   Rcpp::Named("log-normal") = LOGNORMAL,
   Rcpp::Named("single") = SINGLE,

claimsBoot/src/glm.cpp

@@ -20,16 +20,21 @@ int validate_rng_(int n_samples);
 //' @export
 // [[Rcpp::export]]
 Rcpp::NumericVector glmBoot(Rcpp::NumericMatrix triangle, int n_boot,
-                            Rcpp::String boot_type, Rcpp::String opt,
+                            Rcpp::String boot_type, Rcpp::String opt = "null",
                             int seed = 42) {
   triangle = na_to_zero(triangle);
   int n_dev = triangle.cols();
 
   int boot_type_ = key[boot_type];
-  int opt_ = key[opt];
+  int opt_;
+  if (boot_type == "parametric") {
+    opt_ = key[opt];
+  } else {
+    opt_ = 0;
+  }
 
   Rcpp::NumericVector reserve(n_boot);
-  glm_boot(n_dev, triangle.begin(), boot_type_, n_boot, opt_, reserve.begin(),
+  glm_boot(n_dev, triangle.begin(), boot_type_, opt_, n_boot, reserve.begin(),
            seed);
   return (reserve);
 };

claimsBoot/src/mack.cpp

@@ -86,7 +86,8 @@ Rcpp::DataFrame mackSim(Rcpp::NumericMatrix triangle, Rcpp::String sim_type,
     mack_sim(n_dev, triangle.begin(), sim_type_, boot_type_, sim_dist_, n_boot,
              n_conf, m_conf, conf.begin(), fort_res.begin(), show_progress,
              seed);
-    res_list[boot_type_ - 1] = fort_res;
+    int idx = i - boot_types.begin();
+    res_list[idx] = fort_res;
   }
   Rcpp::DataFrame res = mackPost(res_list, boot_types, sim_type);
   return (res);

claimsBoot/src/mod_glm.f90

@@ -17,7 +17,7 @@ module mod_glm
 contains
 
   subroutine sim(n_dev, triangle, sim_type, boot_type_, n_conf, m_conf, conf, &
-      n_boot, res, show_progress, seed) bind(c, name="glm_sim")
+    n_boot, res, show_progress, seed) bind(c, name="glm_sim")
     integer(c_int), intent(in), value :: n_dev, n_boot, n_conf, m_conf
     integer(c_int), intent(in), value :: sim_type, boot_type_
     real(c_double), intent(in) :: conf(n_conf, m_conf)
@@ -72,7 +72,6 @@ subroutine sim(n_dev, triangle, sim_type, boot_type_, n_conf, m_conf, conf, &
     !$omp do schedule(dynamic, 25)
     do i_sim = 1, n_conf
       if (show_progress) call pgbar_incr(pgbar)
-
       if (sim_type == SINGLE) then
         outlier_rowidx = int(conf(i_sim, 1))
         outlier_colidx = int(conf(i_sim, 2))
@@ -176,10 +175,13 @@ subroutine boot(n_dev, triangle, n_boot, reserve)
     integer(c_int) :: n_pts, n_covs, n_pred, n_obs
     real(c_double), allocatable :: betas(:), X_pred(:, :), y_pred(:)
     real(c_double), allocatable :: triangle_boot(:, :), betas_boot(:), resids_boot(:, :)
+    real(c_double), allocatable :: triangle_pred(:, :)
 
     real(c_double) :: disp, disp_boot
     real(c_double) :: lambda, mean, sd, shape, scale
 
+    real(c_double) :: resid_sim
+
     real(c_double), allocatable :: flat_resids(:)
 
     integer(c_int) :: i, j, k, i_boot
@@ -196,6 +198,7 @@ subroutine boot(n_dev, triangle, n_boot, reserve)
 
     allocate(resids_boot(n_pts, n_pts), source = 0._c_double)
     allocate(triangle_boot(n_dev, n_dev), source = 0._c_double)
+    allocate(triangle_pred(n_dev, n_dev), source = 0._c_double)
 
     allocate(X_pred(n_pred, n_covs), source = 0._c_double)
     allocate(y_pred(n_pred), source = 0._c_double)
@@ -222,6 +225,7 @@ subroutine boot(n_dev, triangle, n_boot, reserve)
         end do
 
       else
+        triangle_boot = 0
         do i = 1, n_dev
           do j = 1, n_dev + 1 - i
             if (dist == NORMAL) then
@@ -261,29 +265,57 @@ subroutine boot(n_dev, triangle, n_boot, reserve)
 
       y_pred = exp(matmul(X_pred, betas_boot))
 
-      do i = 1, n_pred
-        lambda = y_pred(i)
-        if (lambda > 1e7) cycle main_loop
-        y_pred(i) = rpois_par(rng, i_thread, lambda)
+      k = 1
+      do i = 2, n_dev
+        do j = n_dev + 1 - i + 1, n_dev
+          triangle_pred(i, j) = y_pred(k)
+          k = k + 1
+        end do
       end do
 
+      if (boot_type == RESID) then
+        k = 1
+        do i = 2, n_dev
+          do j = n_dev + 1 - i + 1, n_dev
+            resid_sim = flat_resids(1 + int(n_pts * runif_par(rng, i_thread)))
+            triangle_pred(i, j) = triangle_pred(i, j) + resid_sim * sqrt(triangle_pred(i, j))
+          end do
+        end do
+      else
+        do i = 1, n_pred
+          lambda = y_pred(i)
+          if (lambda > 1e7) cycle main_loop
+          y_pred(i) = rpois_par(rng, i_thread, lambda)
+        end do
+      end if
+
       reserve(i_boot) = sum(y_pred)
       i_boot = i_boot + 1
     end do main_loop
   end subroutine boot
 
-  subroutine boot_cpp(n_dev, triangle, boot_type_, opt, reserve, n_boot, seed) bind(c, name="glm_boot")
+  subroutine boot_cpp(n_dev, triangle, boot_type_, opt, n_boot, reserve, seed) bind(c, name="glm_boot")
     integer(c_int), intent(in), value :: n_boot, n_dev
     real(c_double), intent(in) :: triangle(n_dev, n_dev)
     integer(c_int), intent(in), value :: boot_type_, opt
     real(c_double), intent(out) :: reserve(n_boot)
     integer(c_int), intent(in), value :: seed
 
+    integer(c_int) :: i, j
+
     boot_type = boot_type_
     if (boot_type_ == PARAM) dist = opt
 
     allocate(mask(n_dev, n_dev))
+    allocate(resids(n_dev, n_dev), source=0._c_double)
+    allocate(triangle_fit(n_dev, n_dev), source=0._c_double)
+
     mask = .true.
+    do j = 1, n_dev
+      do i = n_dev + 2 - j, n_dev
+        mask(i, j) = .false.
+      end do
+    end do
 
     if (first_call) then
       rng = init_rng(1, seed)
@@ -293,6 +325,10 @@ subroutine boot_cpp(n_dev, triangle, boot_type_, opt, reserve, n_boot, seed) bin
     i_thread = 0
 
     call boot(n_dev, triangle, n_boot, reserve)
+
+    deallocate(mask)
+    deallocate(resids)
+    deallocate(triangle_fit)
   end subroutine boot_cpp
 
   subroutine fit(triangle, betas, disp, use_mask, comp, status)
@@ -383,14 +419,6 @@ subroutine fit(triangle, betas, disp, use_mask, comp, status)
       end do
     end do
 
-    call dgetrf(n_obs, n_covs, X, n_obs, IPIV, info)
-    temp = 1
-    do i = 1, n_covs
-      temp = temp * X(i, i)
-    end do
-
-    if (temp == 0) print *, raise(2)
-
     ! Initialise GLM coefficients.
     betas = spread(0._c_double, 1, n_covs)
     mu = y + 0.1
@@ -400,7 +428,7 @@ subroutine fit(triangle, betas, disp, use_mask, comp, status)
     IPIV = 0
     info = 0
     diff = 1E6
-    eps = 1E-3
+    eps = 1E-5
     n_iter = 0
     status = SUCCESS
     do while (diff > eps .and. n_iter < max_iter)
@@ -415,7 +443,10 @@ subroutine fit(triangle, betas, disp, use_mask, comp, status)
       lhs = matmul(W, X)
       rhs = matmul(W, z)
 
-      if (any(isnan(rhs))) print *, raise(2)
+      if (any(isnan(rhs)) .or. any(isnan(lhs))) then
+        status = FAILURE
+        return
+      end if
 
       call dgels('N', n_obs, n_covs, 1, lhs, n_obs, rhs, n_obs, work, lwork, info)
 
@@ -454,8 +485,8 @@ subroutine fit(triangle, betas, disp, use_mask, comp, status)
       k = 1
       do i = 1, n_dev
         do j = 1, n_dev + 1 - i
-            triangle_fit(i, j) = y_fit(k)
-            k = k + 1
+          triangle_fit(i, j) = y_fit(k)
+          k = k + 1
         end do
       end do
 
@@ -562,4 +593,4 @@ end module mod_glm
 
 ! lhs = matmul(matmul(transpose(X), W), X)
 ! rhs = matmul(matmul(transpose(X), W), z)
-! call dgesv(n_covs, 1, lhs, n_covs, IPIV, rhs, n_covs, info)
+! call dgesv(n_covs, 1, lhs, n_covs, IPIV, rhs, n_covs, info)