## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", eval = FALSE # heavy computation — run manually ) ## ----packages----------------------------------------------------------------- # library(rmdcev) # library(apollo) # library(dplyr) # library(tidyr) # library(bench) # for timing; install.packages("bench") if needed ## ----fn-to-wide--------------------------------------------------------------- # to_wide <- function(data_long, nalts) { # data_long <- as.data.frame(data_long) # # # Individual-level covariates (b-cols that don't vary within id) # # b1..b_npsi_j are alt-specific; remaining b-cols are individual-specific # all_b <- grep("^b[0-9]", names(data_long), value = TRUE) # # detect individual-specific cols: same value for all rows of an id # ind_b <- all_b[sapply(all_b, function(v) { # all(tapply(data_long[[v]], data_long$id, function(x) length(unique(x)) == 1)) # })] # # ind_vars <- data_long %>% distinct(id, income, across(all_of(ind_b))) # # alt_vars <- data_long %>% # select(id, alt, quant, price, all_of(setdiff(all_b, ind_b))) %>% # pivot_wider(names_from = alt, values_from = c(quant, price, setdiff(all_b, ind_b))) # # db <- left_join(ind_vars, alt_vars, by = "id") # # # Numeraire quantity # pj_xj <- Reduce("+", lapply(1:nalts, function(j) { # db[[paste0("price_", j)]] * db[[paste0("quant_", j)]] # })) # db$quant_num <- db$income - pj_xj # db # } ## ----fn-apollo-prob----------------------------------------------------------- # make_apollo_prob <- function(nalts, model = c("hybrid0", "gamma"), # npsi_j = 1, ind_col = "b2") { # model <- match.arg(model) # alt_names <- paste0("alt", 1:nalts) # all_names <- c("outside", alt_names) # # function(apollo_beta, apollo_inputs, functionality = "estimate") { # apollo_attach(apollo_beta, apollo_inputs) # on.exit(apollo_detach(apollo_beta, apollo_inputs)) # P <- list() # # # Utility for each non-numeraire alternative # V <- setNames(vector("list", nalts + 1), all_names) # V[["outside"]] <- 0 # for (j in seq_len(nalts)) { # vj <- get(ind_col) * psi_ind # individual-specific component # for (k in seq_len(npsi_j)) { # vj <- vj + get(paste0("b", k, "_", j)) * get(paste0("psi_j", k)) # } # V[[paste0("alt", j)]] <- vj # } # # # Alpha: hybrid0 = 0 for all; gamma = free for outside, 0 for others # alpha_out <- if (model == "gamma") exp(ln_alpha_out) / (1 + exp(ln_alpha_out)) else 0 # alpha <- setNames( # c(list(alpha_out), lapply(seq_len(nalts), function(j) 0)), # all_names # ) # # gamma <- setNames( # c(list(1), lapply(seq_len(nalts), function(j) exp(get(paste0("ln_g", j))))), # all_names # ) # # sigma <- exp(ln_sigma) # # cChoice <- setNames( # c(list(quant_num), lapply(seq_len(nalts), function(j) get(paste0("quant_", j)))), # all_names # ) # cost <- setNames( # c(list(1), lapply(seq_len(nalts), function(j) get(paste0("price_", j)))), # all_names # ) # avail <- setNames( # c(list(1), lapply(seq_len(nalts), function(j) 1)), # all_names # ) # # mdcev_settings <- list( # alternatives = all_names, # budget = income, # V = V, # alpha = alpha, # gamma = gamma, # sigma = sigma, # continuousChoice = cChoice, # cost = cost, # avail = avail # ) # # P[["model"]] <- apollo_mdcev(mdcev_settings, functionality) # P <- apollo_prepareProb(P, apollo_inputs, functionality) # return(P) # } # } ## ----fn-apollo-beta----------------------------------------------------------- # make_apollo_beta <- function(nalts, model = c("hybrid0", "gamma"), npsi_j = 1) { # model <- match.arg(model) # betas <- c(psi_ind = 0) # for (k in seq_len(npsi_j)) betas[paste0("psi_j", k)] <- 0 # for (j in seq_len(nalts)) betas[paste0("ln_g", j)] <- 0 # if (model == "gamma") betas["ln_alpha_out"] <- 0 # betas["ln_sigma"] <- 0 # betas # } ## ----fn-fit-both-------------------------------------------------------------- # fit_both <- function(sim, model = "hybrid0", npsi_j = 1, ind_col = "b2") { # nalts <- sim$data %>% distinct(alt) %>% nrow() # # # ---- rmdcev ---- # rmdcev_formula <- as.formula(paste0( # "~ ", paste(paste0("b", 1:(npsi_j + 1)), collapse = " + ") # )) # t_rmdcev <- system.time( # fit_r <- mdcev( # formula = rmdcev_formula, # data = sim$data, # psi_ascs = 0, # model = model, # algorithm = "MLE", # print_iterations = FALSE, # backend = "rstan" # ) # ) # # # ---- Apollo ---- # database <- to_wide(sim$data, nalts) # db_col <- paste0("b", npsi_j + 1) # individual-specific column # # apollo_control <- list( # modelName = paste0("apollo_", model), # modelDescr = paste0("Apollo ", model), # indivID = "id", # outputDirectory = tempdir(), # nCores = 1 # ) # apollo_beta <- make_apollo_beta(nalts, model, npsi_j) # apollo_fixed <- character(0) # # apollo_inputs <- apollo_validateInputs( # apollo_beta = apollo_beta, # apollo_fixed = apollo_fixed, # database = database, # apollo_control = apollo_control # ) # # apollo_prob <- make_apollo_prob(nalts, model, npsi_j, db_col) # # t_apollo <- system.time( # fit_a <- apollo_estimate( # apollo_beta, apollo_fixed, apollo_prob, apollo_inputs, # estimate_settings = list(writeIter = FALSE, hessianRoutine = "none", # printLevel = 0L) # ) # ) # # # ---- Extract estimates ---- # r_psi <- as.numeric(fit_r$stan_fit$par$psi) # r_gamma <- as.numeric(fit_r$stan_fit$par$gamma) # r_scale <- as.numeric(fit_r$stan_fit$par$scale) # # a_est <- fit_a$estimate # a_gamma <- exp(a_est[grep("^ln_g", names(a_est))]) # a_scale <- exp(a_est["ln_sigma"]) # # # Map column names: psi_j1 = alt-varying, psi_ind = individual-varying # a_psi_j <- a_est[grep("^psi_j", names(a_est))] # a_psi_ind <- a_est["psi_ind"] # a_psi <- c(a_psi_j, a_psi_ind) # # list( # ll = c(rmdcev = fit_r$log.likelihood, apollo = fit_a$LLout), # psi = rbind(rmdcev = r_psi, apollo = a_psi), # gamma = rbind(rmdcev = r_gamma[1:min(5, nalts)], # apollo = a_gamma[1:min(5, nalts)]), # scale = c(rmdcev = r_scale, apollo = as.numeric(a_scale)), # time_s = c(rmdcev = t_rmdcev["elapsed"], apollo = t_apollo["elapsed"]) # ) # } ## ----param-recovery----------------------------------------------------------- # set.seed(42) # true_psi_j <- 0.5 # true_psi_ind <- -1.5 # true_gamma <- rep(5, 5) # true_scale <- 1.0 # # sim_h0 <- GenerateMDCEVData( # model = "hybrid0", # nobs = 500, # nalts = 5, # psi_j_parms = true_psi_j, # psi_i_parms = true_psi_ind, # gamma_parms = true_gamma, # scale_parms = true_scale # ) # # cmp_h0 <- fit_both(sim_h0, model = "hybrid0", npsi_j = 1) # # cat("=== hybrid0 parameter comparison ===\n") # cat(sprintf("LL: rmdcev = %.2f apollo = %.2f\n", # cmp_h0$ll["rmdcev"], cmp_h0$ll["apollo"])) # cat(sprintf(" LL difference = %.6f\n", # abs(cmp_h0$ll["rmdcev"] - cmp_h0$ll["apollo"]))) # # cat("\nPsi (true: psi_j=0.5, psi_ind=-1.5)\n") # print(round(cmp_h0$psi, 4)) # # cat("\nGamma (true: all 5)\n") # print(round(cmp_h0$gamma, 4)) # # cat(sprintf("\nScale (true: 1.0): rmdcev=%.4f apollo=%.4f\n", # cmp_h0$scale["rmdcev"], cmp_h0$scale["apollo"])) ## ----param-recovery-gamma----------------------------------------------------- # set.seed(42) # true_alpha <- 0.5 # # sim_g <- GenerateMDCEVData( # model = "gamma", # nobs = 500, # nalts = 5, # psi_j_parms = true_psi_j, # psi_i_parms = true_psi_ind, # gamma_parms = true_gamma, # alpha_parms = true_alpha, # scale_parms = true_scale # ) # # cmp_g <- fit_both(sim_g, model = "gamma", npsi_j = 1) # # cat("=== gamma parameter comparison ===\n") # cat(sprintf("LL: rmdcev = %.2f apollo = %.2f\n", # cmp_g$ll["rmdcev"], cmp_g$ll["apollo"])) # cat(sprintf(" LL difference = %.6f\n", # abs(cmp_g$ll["rmdcev"] - cmp_g$ll["apollo"]))) # # cat("\nPsi:\n"); print(round(cmp_g$psi, 4)) # cat("\nGamma (first 5):\n"); print(round(cmp_g$gamma, 4)) # cat(sprintf("\nScale: rmdcev=%.4f apollo=%.4f\n", # cmp_g$scale["rmdcev"], cmp_g$scale["apollo"])) ## ----timing-setup------------------------------------------------------------- # grid <- expand.grid( # n = c(500, 1000, 5000), # J = c(50, 100), # rep = 1:3, # stringsAsFactors = FALSE # ) %>% # filter(!(n == 1000 & J == 100)) # keep n=1000 for J=50 only ## ----timing-mle--------------------------------------------------------------- # timing_mle <- lapply(seq_len(nrow(grid)), function(i) { # n <- grid$n[i]; J <- grid$J[i]; r <- grid$rep[i] # set.seed(r * 100 + i) # # sim <- GenerateMDCEVData( # model = "hybrid0", # nobs = n, # nalts = J, # psi_j_parms = 0.5, # psi_i_parms = -1.5, # gamma_parms = rep(5, J), # scale_parms = 1.0 # ) # # res <- fit_both(sim, model = "hybrid0", npsi_j = 1) # data.frame(n=n, J=J, rep=r, # rmdcev_s = res$time_s["rmdcev"], # apollo_s = res$time_s["apollo"], # ll_diff = abs(res$ll["rmdcev"] - res$ll["apollo"])) # }) # # timing_mle <- do.call(rbind, timing_mle) ## ----timing-mle-summary------------------------------------------------------- # summary_mle <- timing_mle %>% # group_by(n, J) %>% # summarise( # rmdcev_mean_s = round(mean(rmdcev_s), 1), # apollo_mean_s = round(mean(apollo_s), 1), # speedup = round(mean(apollo_s) / mean(rmdcev_s), 1), # ll_diff_mean = signif(mean(ll_diff), 3), # .groups = "drop" # ) # # cat("=== MLE estimation timing (seconds, 3 replicates) ===\n") # print(summary_mle, n = Inf) ## ----sim-apollo-prob5--------------------------------------------------------- # # Concrete apollo_probabilities for J=5 (used only in this section) # apollo_probabilities_j5 <- function(apollo_beta, apollo_inputs, functionality = "estimate") { # apollo_attach(apollo_beta, apollo_inputs) # on.exit(apollo_detach(apollo_beta, apollo_inputs)) # P <- list() # # V <- list( # outside = 0, # alt1 = psi_ind * b2 + psi_j1 * b1_1, # alt2 = psi_ind * b2 + psi_j1 * b1_2, # alt3 = psi_ind * b2 + psi_j1 * b1_3, # alt4 = psi_ind * b2 + psi_j1 * b1_4, # alt5 = psi_ind * b2 + psi_j1 * b1_5 # ) # mdcev_settings <- list( # alternatives = c("outside","alt1","alt2","alt3","alt4","alt5"), # budget = income, V = V, # alpha = list(outside=0, alt1=0, alt2=0, alt3=0, alt4=0, alt5=0), # gamma = list(outside=1, alt1=exp(ln_g1), alt2=exp(ln_g2), alt3=exp(ln_g3), # alt4=exp(ln_g4), alt5=exp(ln_g5)), # sigma = exp(ln_sigma), # continuousChoice = list(outside=quant_num, alt1=quant_1, alt2=quant_2, # alt3=quant_3, alt4=quant_4, alt5=quant_5), # cost = list(outside=1, alt1=price_1, alt2=price_2, # alt3=price_3, alt4=price_4, alt5=price_5), # avail = list(outside=1, alt1=1, alt2=1, alt3=1, alt4=1, alt5=1) # ) # P[["model"]] <- apollo_mdcev(mdcev_settings, functionality) # P <- apollo_prepareProb(P, apollo_inputs, functionality) # return(P) # } ## ----sim-convergence---------------------------------------------------------- # set.seed(42) # J5 <- 5 # n500 <- 500 # # sim5 <- GenerateMDCEVData( # model = "hybrid0", # nobs = n500, # nalts = J5, # psi_j_parms = 0.5, # psi_i_parms = -1.5, # gamma_parms = rep(5, J5), # scale_parms = 1.0 # ) # # # ── Fit both packages ────────────────────────────────────────────────────── # fit_r5 <- mdcev( # formula = ~ b1 + b2, data = sim5$data, psi_ascs = 0, # model = "hybrid0", algorithm = "MLE", # print_iterations = FALSE, backend = "rstan" # ) # # database5 <- to_wide(sim5$data, J5) # ac5 <- list(modelName = "apl_sim5", modelDescr = "sim", indivID = "id", # outputDirectory = tempdir(), nCores = 1) # ab5 <- make_apollo_beta(J5, "hybrid0", npsi_j = 1) # ai5 <- suppressMessages( # apollo_validateInputs(apollo_beta = ab5, apollo_fixed = character(0), # database = database5, apollo_control = ac5)) # fit_a5 <- suppressMessages(suppressWarnings( # apollo_estimate(ab5, character(0), apollo_probabilities_j5, ai5, # estimate_settings = list(writeIter = FALSE, hessianRoutine = "none", # printLevel = 0L)) # )) # # # ── 1% price-change policy ───────────────────────────────────────────────── # baseline_prices <- matrix(as.numeric(fit_r5$stan_data$price_j), nrow = n500) # price_change <- 0.01 * colMeans(baseline_prices) # 1% of mean price per alt # # # rmdcev: additive policy vector (0 for numeraire, then per-alt changes) # policies5 <- CreateBlankPolicies(npols = 2, fit_r5, price_change_only = TRUE) # policies5$price_p[[2]] <- c(0, price_change) # # df_sim5 <- PrepareSimulationData(fit_r5, policies5, nsims = 1) # # t_rmdcev5 <- system.time( # demand_r5 <- mdcev.sim( # df_sim5$df_indiv, df_common = df_sim5$df_common, # sim_options = df_sim5$sim_options, # cond_err = 0, nerrs = 100, sim_type = "demand" # ) # ) # # # Extract mean demand per alt: demand[[indiv]][[sim]][policy, alt_index] # mean_demand <- function(demand_list, policy_idx, alt_idx) { # mean(sapply(demand_list, function(ind) # mean(sapply(ind, function(s) s[policy_idx, alt_idx])))) # } # d_base_r <- sapply(2:(J5 + 1), function(a) mean_demand(demand_r5, 1, a)) # d_policy_r <- sapply(2:(J5 + 1), function(a) mean_demand(demand_r5, 2, a)) # delta_r <- d_policy_r - d_base_r # # # Apollo: predict at baseline and policy databases # # Policy database: same additive change as rmdcev (+ 0.01 * mean_price_j per alt) # database5_pol <- database5 # for (j in 1:J5) { # database5_pol[[paste0("price_", j)]] <- database5[[paste0("price_", j)]] + price_change[j] # } # # Recompute quant_num at observed quantities under new prices # pj_xj_pol <- Reduce("+", lapply(1:J5, function(j) # database5_pol[[paste0("price_", j)]] * database5_pol[[paste0("quant_", j)]])) # database5_pol$quant_num <- database5_pol$income - pj_xj_pol # # ai5_pol <- suppressMessages( # apollo_validateInputs(apollo_beta = fit_a5$estimate, apollo_fixed = character(0), # database = database5_pol, apollo_control = ac5)) # # t_apollo5_base <- system.time( # pred_base <- suppressMessages(suppressWarnings( # apollo_prediction(fit_a5, apollo_probabilities_j5, ai5, # prediction_settings = list(nRep = 100, silent = TRUE)) # )) # ) # t_apollo5_pol <- system.time( # pred_pol <- suppressMessages(suppressWarnings( # apollo_prediction(fit_a5, apollo_probabilities_j5, ai5_pol, # prediction_settings = list(nRep = 100, silent = TRUE)) # )) # ) # t_apollo5 <- t_apollo5_base + t_apollo5_pol # # delta_a <- sapply(1:J5, function(j) { # mean(pred_pol[[paste0("alt", j, "_cont_mean")]]) - # mean(pred_base[[paste0("alt", j, "_cont_mean")]]) # }) # # cat("=== Demand change: 1% price increase (n=500, J=5) ===\n") # cat(sprintf("%-12s %8s %8s %8s\n", "Alt", "rmdcev", "Apollo", "Diff")) # for (j in 1:J5) { # cat(sprintf("alt%-9d %8.4f %8.4f %8.4f\n", j, delta_r[j], delta_a[j], # abs(delta_r[j] - delta_a[j]))) # } # cat(sprintf("\nMean abs difference: %.5f\n", mean(abs(delta_r - delta_a)))) # cat(sprintf("rmdcev %% change: %s%%\n", # paste(round(100 * delta_r / d_base_r, 2), collapse = ", "))) # # cat(sprintf("\nTiming: rmdcev %.2f s | Apollo (2 predictions) %.2f s\n", # t_rmdcev5["elapsed"], t_apollo5["elapsed"])) ## ----make-concrete-prob------------------------------------------------------- # # Dynamically generate a concrete apollo_probabilities function for any J. # # apollo_prediction() does source-code inspection and cannot use factory closures, # # so we build and eval() a function body with the right number of alternatives. # make_apollo_prob_concrete <- function(nalts, npsi_j = 1, ind_col = "b2") { # all_str <- paste0('c("outside", ', paste0('"alt', seq_len(nalts), '"', collapse = ", "), ")") # V_items <- paste(paste0("alt", seq_len(nalts), " = psi_ind * ", ind_col, # " + psi_j1 * b1_", seq_len(nalts)), collapse = ",\n ") # alpha_items <- paste(c("outside = 0", paste0("alt", seq_len(nalts), " = 0")), collapse = ", ") # gamma_items <- paste(c("outside = 1", paste0("alt", seq_len(nalts), " = exp(ln_g", seq_len(nalts), ")")), collapse = ", ") # choice_items <- paste(c("outside = quant_num", paste0("alt", seq_len(nalts), " = quant_", seq_len(nalts))), collapse = ", ") # cost_items <- paste(c("outside = 1", paste0("alt", seq_len(nalts), " = price_", seq_len(nalts))), collapse = ", ") # avail_items <- paste(c("outside = 1", paste0("alt", seq_len(nalts), " = 1")), collapse = ", ") # # eval(parse(text = sprintf( # 'function(apollo_beta, apollo_inputs, functionality = "estimate") { # apollo_attach(apollo_beta, apollo_inputs) # on.exit(apollo_detach(apollo_beta, apollo_inputs)) # P <- list() # V <- list(outside = 0, %s) # mdcev_settings <- list( # alternatives = %s, budget = income, V = V, # alpha = list(%s), gamma = list(%s), sigma = exp(ln_sigma), # continuousChoice = list(%s), cost = list(%s), avail = list(%s) # ) # P[["model"]] <- apollo_mdcev(mdcev_settings, functionality) # P <- apollo_prepareProb(P, apollo_inputs, functionality) # return(P) # }', V_items, all_str, alpha_items, gamma_items, choice_items, cost_items, avail_items))) # } ## ----timing-sim-grid---------------------------------------------------------- # timing_sim <- lapply(c(50, 100), function(J) { # n <- 500 # set.seed(J) # # sim <- GenerateMDCEVData( # model = "hybrid0", nobs = n, nalts = J, # psi_j_parms = 0.5, psi_i_parms = -1.5, # gamma_parms = rep(5, J), scale_parms = 1.0 # ) # # # ── rmdcev ────────────────────────────────────────────────────────────── # fit_r <- mdcev(~ b1 + b2, data = sim$data, psi_ascs = 0, model = "hybrid0", # algorithm = "MLE", print_iterations = FALSE, backend = "rstan") # bp <- matrix(as.numeric(fit_r$stan_data$price_j), nrow = n) # pc <- 0.01 * colMeans(bp) # pols <- CreateBlankPolicies(npols = 2, fit_r, price_change_only = TRUE) # pols$price_p[[2]] <- c(0, pc) # df_s <- PrepareSimulationData(fit_r, pols, nsims = 1) # # rmdcev_times <- sapply(1:3, function(r) # system.time(mdcev.sim(df_s$df_indiv, df_common = df_s$df_common, # sim_options = df_s$sim_options, # cond_err = 0, nerrs = 100, sim_type = "demand"))["elapsed"]) # # # ── Apollo ────────────────────────────────────────────────────────────── # database <- to_wide(sim$data, J) # ac <- list(modelName = paste0("apl_", J), modelDescr = "", indivID = "id", # outputDirectory = tempdir(), nCores = 1) # ab <- make_apollo_beta(J, npsi_j = 1) # pfn <- make_apollo_prob_concrete(J) # ai <- apollo_validateInputs(apollo_beta = ab, apollo_fixed = character(0), # database = database, apollo_control = ac) # fit_a <- apollo_estimate(ab, character(0), pfn, ai, # estimate_settings = list(writeIter = FALSE, hessianRoutine = "none", # printLevel = 0L)) # # # Policy database: 1% additive price increase # db_pol <- database # for (j in seq_len(J)) db_pol[[paste0("price_", j)]] <- database[[paste0("price_", j)]] + pc[j] # pj_xj_pol <- Reduce("+", lapply(seq_len(J), function(j) # db_pol[[paste0("price_", j)]] * db_pol[[paste0("quant_", j)]])) # db_pol$quant_num <- db_pol$income - pj_xj_pol # ai_pol <- apollo_validateInputs(apollo_beta = fit_a$estimate, apollo_fixed = character(0), # database = db_pol, apollo_control = ac) # # apollo_times <- sapply(1:3, function(r) { # t1 <- system.time(apollo_prediction(fit_a, pfn, ai, # prediction_settings = list(nRep = 100, silent = TRUE)))["elapsed"] # t2 <- system.time(apollo_prediction(fit_a, pfn, ai_pol, # prediction_settings = list(nRep = 100, silent = TRUE)))["elapsed"] # t1 + t2 # }) # # data.frame(n = n, J = J, # rmdcev_mean = round(mean(rmdcev_times), 2), # rmdcev_sd = round(sd(rmdcev_times), 2), # apollo_mean = round(mean(apollo_times), 2), # apollo_sd = round(sd(apollo_times), 2), # speedup = round(mean(apollo_times) / mean(rmdcev_times), 1)) # }) # # timing_sim <- do.call(rbind, timing_sim) # print(timing_sim, row.names = FALSE) ## ----session------------------------------------------------------------------ # sessionInfo()