This R Markdown document contains the code for performing simulation studies related to the present method for setting acceptance criteria.
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.8 ✔ dplyr 1.0.9
## ✔ tidyr 1.2.0 ✔ stringr 1.4.0
## ✔ readr 2.1.2 ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(furrr)
## Loading required package: future
##
## Attaching package: 'future'
##
## The following object is masked from 'package:rmarkdown':
##
## run
library(cmstatr)
library(tikzDevice)
plan(multisession, workers = 12)
factors <- read_csv(
"factors.csv",
col_types = cols(
alpha = col_double(),
n = col_double(),
m = col_double(),
k1 = col_double(),
k2 = col_double()
))
mu_pop <- 0
sd_pop <- 1
n_qual_1 <- 5000
n_equiv_1 <- 5000
n_qual_2 <- 2500
n_equiv_2 <- 2500
sim_alpha <- 0.05
delta_mean <- seq(0, 2, length.out = 9)
delta_sd <- seq(1, 5, length.out = 9)
delta_min_indiv <- seq(0, 4, length.out = 9)
delta_mixture <- seq(0, 4, length.out = 9)
pi.mixture <- 0.4
sim_n_values <- unique(factors$n[factors$n <= 100])
sim_m_values <- unique(factors$m[factors$m %% 2 == 0])
k_factors <- unique(factors$m) %>%
map_dfr(function(m) {
k <- k_equiv(sim_alpha, m)
tibble(
m = m,
k1 = k[1],
k2 = k[2]
)
})
sim_equiv <- map_dfr(
sim_m_values,
function(cur_m) {
k1 <- k_factors$k1[k_factors$m == cur_m]
k2 <- k_factors$k2[k_factors$m == cur_m]
future_map_dfr(
sim_n_values,
.options = furrr_options(seed = 1234),
function(cur_n) {
r1 <- factors$k1[factors$n == cur_n & factors$m == cur_m
& factors$alpha == sim_alpha]
r2 <- factors$k2[factors$n == cur_n & factors$m == cur_m
& factors$alpha == sim_alpha]
equiv_samples <- future_map_dfr(
1:n_equiv_1,
.options = furrr_options(seed = 2345),
function(i_equiv) {
x <- rnorm(cur_m, mu_pop, sd_pop)
data.frame(
min = min(x),
avg = mean(x)
)
}
)
future_map_dfr(
1:n_qual_1,
.options = furrr_options(seed = 3456),
function(i_qual) {
x_qual <- rnorm(cur_n, mu_pop, sd_pop)
avg_qual <- mean(x_qual)
sd_qual <- sd(x_qual)
rbind(
data.frame(
method = factor("Vangel", c("Two-Sample", "Vangel")),
accept_min = equiv_samples$min > avg_qual - k1 * sd_qual,
accept_avg = equiv_samples$avg > avg_qual - k2 * sd_qual
),
data.frame(
method = factor("Two-Sample", c("Two-Sample", "Vangel")),
accept_min = equiv_samples$min > avg_qual - r1 * sd_qual,
accept_avg = equiv_samples$avg > avg_qual - r2 * sd_qual
)
)
}
) %>%
mutate(accept = accept_min & accept_avg) %>%
group_by(method) %>%
summarise(`Rejection Rate` = sum(!accept) / n(),
.groups = "drop") %>%
mutate(n = cur_n, m = cur_m)
}
)
}
)
# equiv_deviation has the arguments:
# i_equiv, delta, cur_m
power_simulation <- function(delta_list, qual_generation, equiv_generation) {
map_dfr(
sim_m_values,
function(cur_m) {
k1 <- k_factors$k1[k_factors$m == cur_m]
k2 <- k_factors$k2[k_factors$m == cur_m]
future_map_dfr(
sim_n_values,
.options = furrr_options(seed = 12345),
function(cur_n) {
r1 <- factors$k1[factors$n == cur_n & factors$m == cur_m
& factors$alpha == sim_alpha]
r2 <- factors$k2[factors$n == cur_n & factors$m == cur_m
& factors$alpha == sim_alpha]
t_crit <- qt(1 - sim_alpha / 2, cur_n - 1)
f_crit <- qf(1 - sim_alpha / 2, cur_m - 1, cur_n - 1)
future_map_dfr(
delta_list,
.options = furrr_options(seed = 23456),
function(delta) {
equiv_samples <- future_map_dfr(
1:n_equiv_2,
.options = furrr_options(seed = 34567),
function(i_equiv) {
x <- equiv_generation(i_equiv, delta, cur_m)
data.frame(
min = min(x),
avg = mean(x),
sd_eq = sd(x)
)
}
)
future_map_dfr(
1:n_qual_2,
.options = furrr_options(seed = 45678),
function(i_qual) {
x_qual <- qual_generation(cur_n)
avg_qual <- mean(x_qual)
sd_qual <- sd(x_qual)
# Ref. Escobar, pp. 55
t_min <- avg_qual - t_crit * sqrt(1 / cur_n + 1 / cur_m) * sd_qual
# Ref. Escobar, pp. 58
f_max <- sd_qual * sqrt(f_crit)
accept_min <- equiv_samples$min > avg_qual - r1 * sd_qual
accept_avg <- equiv_samples$avg > avg_qual - r2 * sd_qual
accept_t <- equiv_samples$avg >= t_min
accept_f <- equiv_samples$sd_eq <= f_max
accept_vangel_min <- equiv_samples$min > avg_qual - k1 * sd_qual
accept_vangel_avg <- equiv_samples$avg > avg_qual - k2 * sd_qual
data.frame(
accept = c(accept_min & accept_avg,
accept_t & accept_f,
accept_vangel_min & accept_vangel_avg),
method = as.factor(c(rep_len("Two-Sample", n_equiv_2),
rep_len("MSD", n_equiv_2),
rep_len("Vangel", n_equiv_2)))
)
}
) %>%
group_by(method) %>%
summarise(`Rejection Rate` = sum(!accept) / n(),
.groups = "drop") %>%
mutate(n = cur_n, m = cur_m, delta = delta)
}
)
}
)
}
)
}
sim_power_mean <- power_simulation(
delta_mean,
function(cur_n) {
rnorm(cur_n, mu_pop, sd_pop)
},
function(i_equiv, delta, cur_m) {
rnorm(cur_m, mu_pop - delta * sd_pop, sd_pop)
}
)
sim_power_sd <- power_simulation(
delta_sd,
function(cur_n) {
rnorm(cur_n, mu_pop, sd_pop)
},
function(i_equiv, delta, cur_m) {
rnorm(cur_m, mu_pop, sd_pop * delta)
}
)
sim_power_min_indiv <- power_simulation(
delta_min_indiv,
function(cur_n) {
rnorm(cur_n, mu_pop, sd_pop)
},
function(i_equiv, delta, cur_m) {
c(
rnorm(cur_m - 1, mu_pop, sd_pop),
rnorm(1, mu_pop - delta * sd_pop, sd_pop)
)
}
)
sim_power_mixture <- power_simulation(
delta_mixture,
function(cur_n) {
rnorm(cur_n, mu_pop, sd_pop)
},
function(i_equiv, delta, cur_m) {
m.1 <- rbinom(1, cur_m, pi.mixture)
c(
rnorm(cur_m - m.1, mu_pop, sd_pop),
rnorm(m.1, mu_pop - delta * sd_pop, sd_pop)
)
}
)
For a lognormal distribution, the coefficient of variation (CV) is given by:
\[ CV = \sqrt{\exp\left(\sigma^2\right) - 1} \]
This can be rearranged to solve for \(\sigma\) as follows. We’ll be setting a CV that is typical of composite materials:
\[ \sigma = \sqrt{\log\left(CV^2 + 1\right)} \]
The mean for a lognormal distribution is given by:
\[ mean = m \exp\left(\frac{1}{2}\sigma^2\right) \]
We’ll set the mean equal to 100.
The standard deviation is given by:
\[ SD = m \exp\left(\frac{1}{2}\sigma^2\right)\sqrt{\exp\left(\sigma^2\right)-1} \]
We’ll set the acceptance mean so that it’s \(\delta\) times the SD lower than the mean used for the qualification data.
\[ mean_{accept} = mean_{pop} - \delta SD \\ m_{accept} \exp\left(\frac{1}{2}\sigma^2\right) = m_{pop} \exp\left(\frac{1}{2}\sigma^2\right) - \delta m_{pop} \exp\left(\frac{1}{2}\sigma^2\right)\sqrt{\exp\left(\sigma^2\right)-1} \\ m_{accept} = m_{pop} - \delta m_{pop} \sqrt{\exp\left(\sigma^2\right)-1} \]
sigma_pop <- sqrt(log(0.06 ^ 2 + 1))
m_pop <- 100 / exp(1/2 * sigma_pop ^ 2)
delta_lognorm <- seq(0, 2, length.out = 9)
The PDFs of the acceptance samples in this simulation for \(\delta=0\), \(\delta=1\) and \(\delta=2\)
ggplot(data.frame(x = c(80 , 120)), aes(x = x)) +
stat_function(fun = dlnorm, args = list(
meanlog = log(m_pop),
sdlog = sigma_pop)) +
stat_function(fun = dlnorm, args = list(
meanlog = log(m_pop - 1.0 * m_pop * sqrt(exp(sigma_pop^2) - 1)),
sdlog = sigma_pop)) +
stat_function(fun = dlnorm, args = list(
meanlog = log(m_pop - 2.0 * m_pop * sqrt(exp(sigma_pop^2) - 1)),
sdlog = sigma_pop))
sim_power_lognormal <- power_simulation(
delta_lognorm,
function(cur_n) {
rlnorm(cur_n, log(m_pop), sigma_pop)
},
function(i_equiv, delta, cur_m) {
rlnorm(cur_n, log(m_pop - delta * m_pop * sqrt(exp(sigma_pop^2) - 1)), sigma_pop)
}
)
For a Weibull distribution with scale parameter \(\eta\) and shape parameter \(\beta\), the CV is given by:
\[ CV = \left\{ \frac{ \Gamma\left[\left(\beta + 2\right) / \beta\right] }{ \left\{\Gamma\left[\left(\beta + 1\right) / \beta\right]\right\}^2 } -1 \right\}^\frac{1}{2} \]
shape_pop <- uniroot(
function(b) {
cv <- 0.06
sqrt(
gamma((b + 2) / b) /
(gamma((b + 1) / b)) ^ 2
- 1
) - cv
},
c(0.1, 1000)
)$root
The mean is given by:
\[ mean = \eta \Gamma\left[\left(\beta + 1\right) / \beta\right] \]
scale_pop <- 100 / gamma((shape_pop + 1) / shape_pop)
The standard deviation for a Weibull distribution is given by:
\[ SD = \eta \sqrt{ \Gamma\left[\left(\beta + 2\right) / \beta\right] - \left\{\Gamma\left[\left(\beta + 1\right) / \beta\right]\right\}^2 } \]
The mean of the acceptance sample is set as follows:
\[ mean_{accept} = mean_{pop} - \delta SD \\ \eta_{accept} \Gamma\left[\left(\beta + 1\right) / \beta\right] = \eta_{pop} \Gamma\left[\left(\beta + 1\right) / \beta\right] - \delta \eta_{pop} \sqrt{ \Gamma\left[\left(\beta + 2\right) / \beta\right] - \left\{\Gamma\left[\left(\beta + 1\right) / \beta\right]\right\}^2 } \\ \eta_{accept} = \eta_{pop} \left(1 - \delta \frac{ \sqrt{ \Gamma\left[\left(\beta + 2\right) / \beta\right] - \left\{\Gamma\left[\left(\beta + 1\right) / \beta\right]\right\}^2 } } { \Gamma\left[\left(\beta + 1\right) / \beta\right] } \right) \]
delta_weibull <- seq(0, 2, length.out = 9)
The PDFs of the acceptance samples in this simulation for \(\delta=0\), \(\delta=1\) and \(\delta=2\)
ggplot(data.frame(x = c(80 , 120)), aes(x = x)) +
stat_function(fun = dweibull, args = list(
shape = shape_pop,
scale = scale_pop)) +
stat_function(fun = dweibull, args = list(
shape = shape_pop,
scale = scale_pop * (1 - 1.0 *(
sqrt(
gamma((shape_pop + 2) / shape_pop) -
(gamma((shape_pop + 1) / shape_pop)) ^ 2
) / gamma((shape_pop + 1) / shape_pop)
)))) +
stat_function(fun = dweibull, args = list(
shape = shape_pop,
scale = scale_pop * (1 - 2.0 *(
sqrt(
gamma((shape_pop + 2) / shape_pop) -
(gamma((shape_pop + 1) / shape_pop)) ^ 2
) / gamma((shape_pop + 1) / shape_pop)
))))
sim_power_weibull <- power_simulation(
delta_weibull,
function(cur_n) {
rweibull(cur_n, shape = shape_pop, scale = scale_pop)
},
function(i_equiv, delta, cur_m) {
rweibull(
cur_n,
shape = shape_pop,
scale = scale_pop * (1 - delta *(
sqrt(
gamma((shape_pop + 2) / shape_pop) -
(gamma((shape_pop + 1) / shape_pop)) ^ 2
) / gamma((shape_pop + 1) / shape_pop)
)))
}
)
save(
mu_pop,
sd_pop,
n_qual_1,
n_equiv_1,
n_qual_2,
n_equiv_2,
sim_alpha,
delta_mean,
delta_sd,
delta_min_indiv,
sim_equiv,
sim_power_mean,
sim_power_sd,
sim_power_min_indiv,
sim_power_mixture,
sim_power_lognormal,
sim_power_weibull,
file = "sim_results.RData"
)
sessionInfo()
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
##
## Random number generation:
## RNG: L'Ecuyer-CMRG
## Normal: Inversion
## Sample: Rejection
##
## locale:
## [1] LC_CTYPE=en_CA.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_CA.UTF-8 LC_COLLATE=en_CA.UTF-8
## [5] LC_MONETARY=en_CA.UTF-8 LC_MESSAGES=en_CA.UTF-8
## [7] LC_PAPER=en_CA.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] tikzDevice_0.12.3.1 cmstatr_0.9.1 furrr_0.3.0
## [4] future_1.27.0 forcats_0.5.1 stringr_1.4.0
## [7] dplyr_1.0.9 purrr_0.3.4 readr_2.1.2
## [10] tidyr_1.2.0 tibble_3.1.8 ggplot2_3.3.6
## [13] tidyverse_1.3.2 rmarkdown_2.14
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.3 sass_0.4.2 bit64_4.0.5
## [4] vroom_1.5.7 jsonlite_1.8.0 modelr_0.1.8
## [7] bslib_0.4.0 assertthat_0.2.1 highr_0.9
## [10] googlesheets4_1.0.0 cellranger_1.1.0 yaml_2.3.5
## [13] globals_0.15.1 pillar_1.8.0 backports_1.4.1
## [16] glue_1.6.2 digest_0.6.29 rvest_1.0.2
## [19] colorspace_2.0-3 htmltools_0.5.3 pkgconfig_2.0.3
## [22] broom_1.0.0 listenv_0.8.0 haven_2.5.0
## [25] scales_1.2.0 tzdb_0.3.0 googledrive_2.0.0
## [28] generics_0.1.3 farver_2.1.1 ellipsis_0.3.2
## [31] cachem_1.0.6 withr_2.5.0 cli_3.3.0
## [34] magrittr_2.0.3 crayon_1.5.1 readxl_1.4.0
## [37] evaluate_0.15 fs_1.5.2 fansi_1.0.3
## [40] parallelly_1.32.1 MASS_7.3-58 xml2_1.3.3
## [43] SuppDists_1.1-9.7 tools_4.2.1 hms_1.1.1
## [46] gargle_1.2.0 lifecycle_1.0.1 munsell_0.5.0
## [49] reprex_2.0.1 kSamples_1.2-9 compiler_4.2.1
## [52] jquerylib_0.1.4 rlang_1.0.4 grid_4.2.1
## [55] filehash_2.4-3 labeling_0.4.2 gtable_0.3.0
## [58] codetools_0.2-18 DBI_1.1.3 R6_2.5.1
## [61] lubridate_1.8.0 knitr_1.39 fastmap_1.1.0
## [64] bit_4.0.4 utf8_1.2.2 stringi_1.7.8
## [67] parallel_4.2.1 vctrs_0.4.1 dbplyr_2.2.1
## [70] tidyselect_1.1.2 xfun_0.31