Added MLE and results tables

Rucknium · Dec 8, 2021 · 54c62a9 · 54c62a9
1 parent e8ac451
commit 54c62a9
Show file tree

Hide file tree

Showing 4 changed files with 89 additions and 1 deletion.
diff --git a/R/OSPEAD-dry-run.R b/R/OSPEAD-dry-run.R
@@ -15,20 +15,32 @@ theta_i <- merge(data.frame(Var1 = 1:max(theta_i)), theta_i, all.x = TRUE)
 theta_i <- theta_i$Freq
 theta_i[is.na(theta_i)] <- 0
 support <- 1:max(xmr.Moser$spend_times_blocks)
+spend_times_blocks <- log(xmr.Moser$spend_times_blocks + 1)
+# offset is 2 to avoid zeroes later on
 
 
+param.trans <- list()
+
 f_D.gamma <- function(param, support) {
   stats::dgamma(log(support), shape = exp(param[1]), rate = exp(param[2]))
 }
 
+param.trans$gamma <- c(exp, exp)
+
 f_D.lnorm <- function(param, support) {
   stats::dlnorm(log(support), meanlog = param[1], sdlog = exp(param[2]))
 }
 
+param.trans$lnorm <- c(I, exp)
+
 f_D.f <- function(param, support) {
   stats::df(log(support), df1 = exp(param[1]), df2 = exp(param[2]), ncp = exp(param[3]) )
 }
 
+param.trans$f <- c(exp, exp, exp)
+
+
+
 L_FGT <- function(param, f_D, support, flavor = 1) {
   alpha <- flavor
   a_i <- f_D(param, support)
@@ -54,16 +66,30 @@ L_Welfare <- function(param, f_D, support, flavor = 1) {
     , na.rm = TRUE )
 }
 
+L_MLE <- function(param, f_D, ...) {
+  (-1) * sum(
+  log( f_D(param, spend_times_blocks + 1) )
+    , na.rm = TRUE)
+}
+
+
+
+
 
 run.iters <- expand.grid(
   f_D = c(f_D.gamma = f_D.gamma, f_D.lnorm = f_D.lnorm, f_D.f = f_D.f), 
   flavor = 1:2,
-  L = c(L_FGT = L_FGT, L_Welfare = L_Welfare)
+  L = c(L_FGT = L_FGT, L_Welfare = L_Welfare, L_MLE = L_MLE)
 )
 
 run.iters$flavor[names(run.iters$L) == "L_Welfare"] <- 
   rep(c(1, 5), each = sum(names(run.iters$L) == "L_Welfare")/2)
 
+run.iters$flavor[names(run.iters$L) == "L_MLE"] <- 0
+# MLE has only one flavor
+
+run.iters <- unique(run.iters)
+
 start.params <- list(
   f_D.gamma = c(1, 1), 
   f_D.lnorm = c(0, 2), 
@@ -108,6 +134,43 @@ for (i in 1:nrow(run.iters.results)) {
   ]]$value
 }
 
+families <- colnames(run.iters.results)[grepl("f_D", colnames(run.iters.results))]
+
+for (family in families) {
+  run.iters.results[run.iters.results$L == "L_MLE", family] <- 
+    2 * length(param.trans[[ gsub("f_D[.]", "", family) ]]) +
+    2 * run.iters.results[run.iters.results$L == "L_MLE", family]
+}
+# AIC formula. Note that the value is the neagtive of the 
+# log likelihood already, so just need to add
+
+
+
+write.csv(run.iters.results, "tables/dry-run/performance.csv", row.names = FALSE)
+
+
+
+minimizer.params <- data.frame(f_D = names(run.iters[[1]]), L = names(run.iters[[3]]), flavor = run.iters[[2]],
+  param_1 = NA, param_2 = NA, param_3 = NA, stringsAsFactors = FALSE)
+
+minimizer.params$f_D <- gsub("f_D[.]", "", minimizer.params$f_D)
+
+for ( i in 1:nrow(minimizer.params)) {
+
+  temp.param.value <- results[[i]]$par
+
+  for ( j in 1:length(temp.param.value)) {
+    temp.param.value[j] <- param.trans[[minimizer.params[i, "f_D"]]][[j]](temp.param.value[j])
+  }
+
+  minimizer.params[i, paste0("param_", seq_along(results[[i]]$par))] <- temp.param.value
+
+}
+
+minimizer.params <- minimizer.params[order(minimizer.params$f_D), ]
+
+write.csv(minimizer.params, "tables/dry-run/minimizer-params.csv", row.names = FALSE)
+
 
 run.iters.obj.fn <- unique(run.iters[, c("L", "flavor")])
 run.iters.obj.fn$L <- names(run.iters.obj.fn$L)
@@ -123,6 +186,9 @@ run.iters.obj.fn$title[run.iters.obj.fn$L == "L_Welfare" & run.iters.obj.fn$flav
 run.iters.obj.fn$title[run.iters.obj.fn$L == "L_Welfare" & run.iters.obj.fn$flavor == 5][[1]] <-
   bquote(L[Welfare][eta] ~ "," ~ eta == 5)
 
+run.iters.obj.fn$title[run.iters.obj.fn$L == "L_MLE"][[1]] <-
+  bquote(L[MLE])
+
 
 distn.name.converter <- c(
   "f_D.gamma" = "Gamma",

diff --git a/images/dry-run/L_MLE-flavor-0.png b/images/dry-run/L_MLE-flavor-0.png
diff --git a/tables/dry-run/minimizer-params.csv b/tables/dry-run/minimizer-params.csv
@@ -0,0 +1,16 @@
+"f_D","L","flavor","param_1","param_2","param_3"
+"f","L_FGT",1,7.73308554515871,118.063330027606,18.2981654946987
+"f","L_FGT",2,11.5921594797106,53.2853002761615,27.3406960709529
+"f","L_Welfare",1,11.9743452905718,59.1764924789675,29.717104784749
+"f","L_Welfare",5,12.8427809452998,33.244205814352,25.8240555862385
+"f","L_MLE",0,48.4622052345434,3443646224.48591,51.0433803223784
+"gamma","L_FGT",1,6.32941654097075,1.87163172217098,NA
+"gamma","L_FGT",2,5.06661209503492,1.62582071865596,NA
+"gamma","L_Welfare",1,5.79170334040451,1.76046364010705,NA
+"gamma","L_Welfare",5,4.89150486745901,1.28402541668774,NA
+"gamma","L_MLE",0,32.4413857250808,15.8015684774846,NA
+"lnorm","L_FGT",1,1.32804852566472,0.295921100314083,NA
+"lnorm","L_FGT",2,1.23371839168609,0.324481375069784,NA
+"lnorm","L_Welfare",1,1.27141411713503,0.311455783458507,NA
+"lnorm","L_Welfare",5,0.5250732421875,0.511099288016679,NA
+"lnorm","L_MLE",0,0.703854078358836,0.18061269184726,NA
diff --git a/tables/dry-run/performance.csv b/tables/dry-run/performance.csv
@@ -0,0 +1,6 @@
+"L","flavor","f_D.gamma","f_D.lnorm","f_D.f"
+"L_FGT",1,0.114593035422099,0.118442304051464,0.112173839576946
+"L_FGT",2,0.0571138439873846,0.0497370994664305,0.0474519635595684
+"L_Welfare",1,0.195514126489678,0.17438112408889,0.166604226378729
+"L_Welfare",5,38.335672781811,9.80396337496936,2.43406002776894
+"L_MLE",0,3072672.94314182,3255804.14458042,3036829.65623237