|
1 | | -#' Gold section search algorithm |
2 | | -#' |
3 | | -#' Implementation of univariate gold-section search algorithm. Algorithm |
4 | | -#' assumes that the objective function has one unique minimum (unimodal) |
5 | | -#' within the specified search interval. Supports integer searches for input |
6 | | -#' parameter. |
7 | | -#' |
8 | | -#' @param f function to be minimized, where the first argument corresponds |
9 | | -#' to the parameter of interest |
10 | | -#' |
11 | | -#' @param interval interval with which to search, specified as |
12 | | -#' \code{c(lower, upper)} |
13 | | -#' |
14 | | -#' @param integer logical; should the parameters be constrained as |
15 | | -#' integer values? |
16 | | -#' |
17 | | -#' @param tol convergence tolerance. Ignored when \code{integer = TRUE} |
18 | | -#' (converges when lowest integer is found between competing boundaries) |
19 | | -#' |
20 | | -#' @param maxiter maximum number of iterations |
21 | | -#' |
22 | | -#' @param ... additional arguments passed to \code{f} |
23 | | -#' |
24 | | -#' @return a \code{list} of the converge and solution information |
25 | | -#' |
26 | | -#' @references |
27 | | -#' |
28 | | -#' Chalmers, R. P., & Adkins, M. C. (2020). Writing Effective and Reliable Monte Carlo Simulations |
29 | | -#' with the SimDesign Package. \code{The Quantitative Methods for Psychology, 16}(4), 248-280. |
30 | | -#' \doi{10.20982/tqmp.16.4.p248} |
31 | | -#' |
32 | | -#' @export |
33 | | -#' |
34 | | -#' @author Phil Chalmers \email{rphilip.chalmers@@gmail.com} |
35 | | -#' @examples |
36 | | -#' |
37 | | -#' y <- c(10, 20, 50) |
38 | | -#' y <- y - mean(y) |
39 | | -#' x <- 1:3 |
40 | | -#' f <- function(beta, y, x) sum( (y - beta*x)^2) |
41 | | -#' |
42 | | -#' betas <- seq(0, 10, length.out=1000) |
43 | | -#' fs <- sapply(betas, f, y=y, x=x) |
44 | | -#' plot(betas, fs, type = 'l', las=1) |
45 | | -#' |
46 | | -#' res <- gold_section(f, interval = c(0,10), y=y, x=x) |
47 | | -#' res |
48 | | -#' |
49 | | -#' # if beta must be an integer |
50 | | -#' betas <- -10:20 |
51 | | -#' fs <- sapply(betas, f, y=y, x=x) |
52 | | -#' plot(betas, fs, las=1) |
53 | | -#' |
54 | | -#' res <- gold_section(f, interval = c(-10,20), integer=TRUE, y=y, x=x) |
55 | | -#' res |
56 | | -#' |
57 | | -#' |
58 | | -#' ############################################# |
59 | | -#' # ANOVA example with Cohen's f= .2, incorporating cost function using SimSolve() |
60 | | -#' # |
61 | | -#' # Goal: solve N (within) and N.groups (between) given the cost function |
62 | | -#' # cost = 5 * N.g + 50 * N.groups |
63 | | -#' # |
64 | | -#' # Setup allows N to be solved with SimSolve() given N.groups and target, |
65 | | -#' # golden_section() used to solve outer N.groups component given cost function |
66 | | -#' |
67 | | -#' # Design row conditions must have N and N.groups |
68 | | -#' costfun <- function(condition) |
69 | | -#' with(condition, 5 * N * N.groups + 50 * N.groups) |
70 | | -#' |
71 | | -#' Generate <- function(condition, fixed_objects = NULL){ |
72 | | -#' Attach(condition) |
73 | | -#' group <- rep(letters[1:N.groups], each=N) |
74 | | -#' means <- rnorm(N.groups, sd=.2) |
75 | | -#' DV <- as.vector(sapply(1:N.groups, \(i) rnorm(N, mean=means[i]))) |
76 | | -#' dat <- data.frame(group, DV) |
77 | | -#' dat |
78 | | -#' } |
79 | | -#' |
80 | | -#' Analyse <- function(condition, dat, fixed_objects = NULL) { |
81 | | -#' mod <- aov(DV ~ group, dat) |
82 | | -#' p <- summary(mod)[[1]]$`Pr(>F)`[1L] |
83 | | -#' p |
84 | | -#' } |
85 | | -#' |
86 | | -#' Summarise <- function(condition, results, fixed_objects = NULL) { |
87 | | -#' ret <- EDR(results) |
88 | | -#' ret |
89 | | -#' } |
90 | | -#' |
91 | | -#' \dontrun{ |
92 | | -#' |
93 | | -#' # Solutions given fixed N.groups input |
94 | | -#' Design <- createDesign(N = NA, N.groups=5:25) |
95 | | -#' |
96 | | -#' # test cost function |
97 | | -#' condition <- Design[1,] |
98 | | -#' condition$N <- 200 |
99 | | -#' costfun(condition) |
100 | | -#' |
101 | | -#' # solve for 1-beta = .95 power |
102 | | -#' results <- SimSolve(design=Design, b=.95, interval=c(10, 500), |
103 | | -#' generate=Generate, analyse=Analyse, summarise=Summarise) |
104 | | -#' results |
105 | | -#' |
106 | | -#' # (optional) confirm results accuracy |
107 | | -#' if(FALSE){ |
108 | | -#' cres <- results |
109 | | -#' cres$N <- ceiling(cres$N) |
110 | | -#' confirm <- runSimulation(cres, replications = 10000, |
111 | | -#' generate=Generate, parallel=TRUE, |
112 | | -#' analyse=Analyse, summarise=Summarise) |
113 | | -#' confirm |
114 | | -#' } |
115 | | -#' |
116 | | -#' costs <- numeric(nrow(results)) |
117 | | -#' for(i in 1:length(costs)) |
118 | | -#' costs[i] <- costfun(ceiling(results[i,])) |
119 | | -#' results_costs <- cbind(results, costs) |
120 | | -#' results_costs |
121 | | -#' |
122 | | -#' library(ggplot2) |
123 | | -#' ggplot(results_costs, aes(N, costs, size = N.groups)) + |
124 | | -#' geom_point() |
125 | | -#' |
126 | | -#' |
127 | | -#' # solved via gold_section() with integer option for N.groups |
128 | | -#' Design <- createDesign(N = NA, N.groups = NA) |
129 | | -#' |
130 | | -#' fn <- function(par, condition){ |
131 | | -#' condition$N.groups <- par |
132 | | -#' res <- SimSolve(design=condition, b=.95, interval=c(10, 500), |
133 | | -#' generate=Generate, analyse=Analyse, summarise=Summarise) |
134 | | -#' costfun(res) |
135 | | -#' } |
136 | | -#' |
137 | | -#' res <- gold_section(fn, interval=c(5,25), integer=TRUE, condition=Design) |
138 | | -#' res |
139 | | -#' |
140 | | -#' } |
| 1 | +# Gold section search algorithm |
| 2 | +# |
| 3 | +# Implementation of univariate gold-section search algorithm. Algorithm |
| 4 | +# assumes that the objective function has one unique minimum (unimodal) |
| 5 | +# within the specified search interval. Supports integer searches for input |
| 6 | +# parameter. |
| 7 | +# |
| 8 | +# @param f function to be minimized, where the first argument corresponds |
| 9 | +# to the parameter of interest |
| 10 | +# |
| 11 | +# @param interval interval with which to search, specified as |
| 12 | +# \code{c(lower, upper)} |
| 13 | +# |
| 14 | +# @param integer logical; should the parameters be constrained as |
| 15 | +# integer values? |
| 16 | +# |
| 17 | +# @param tol convergence tolerance. Ignored when \code{integer = TRUE} |
| 18 | +# (converges when lowest integer is found between competing boundaries) |
| 19 | +# |
| 20 | +# @param maxiter maximum number of iterations |
| 21 | +# |
| 22 | +# @param ... additional arguments passed to \code{f} |
| 23 | +# |
| 24 | +# @return a \code{list} of the converge and solution information |
| 25 | +# |
| 26 | +# @references |
| 27 | +# |
| 28 | +# Chalmers, R. P., & Adkins, M. C. (2020). Writing Effective and Reliable Monte Carlo Simulations |
| 29 | +# with the SimDesign Package. \code{The Quantitative Methods for Psychology, 16}(4), 248-280. |
| 30 | +# \doi{10.20982/tqmp.16.4.p248} |
| 31 | +# |
| 32 | +# @export |
| 33 | +# |
| 34 | +# @author Phil Chalmers \email{rphilip.chalmers@@gmail.com} |
| 35 | +# @examples |
| 36 | +# |
| 37 | +# y <- c(10, 20, 50) |
| 38 | +# y <- y - mean(y) |
| 39 | +# x <- 1:3 |
| 40 | +# f <- function(beta, y, x) sum( (y - beta*x)^2) |
| 41 | +# |
| 42 | +# betas <- seq(0, 10, length.out=1000) |
| 43 | +# fs <- sapply(betas, f, y=y, x=x) |
| 44 | +# plot(betas, fs, type = 'l', las=1) |
| 45 | +# |
| 46 | +# res <- gold_section(f, interval = c(0,10), y=y, x=x) |
| 47 | +# res |
| 48 | +# |
| 49 | +# # if beta must be an integer |
| 50 | +# betas <- -10:20 |
| 51 | +# fs <- sapply(betas, f, y=y, x=x) |
| 52 | +# plot(betas, fs, las=1) |
| 53 | +# min(fs) |
| 54 | +# |
| 55 | +# res <- gold_section(f, interval = c(-10,20), integer=TRUE, y=y, x=x) |
| 56 | +# res |
| 57 | +# |
| 58 | +# |
| 59 | +# ############################################# |
| 60 | +# # ANOVA example with Cohen's f= .2, incorporating cost function using SimSolve() |
| 61 | +# # |
| 62 | +# # Goal: solve N (within) and N.groups (between) given the cost function |
| 63 | +# # cost = 5 * N.g + 50 * N.groups |
| 64 | +# # |
| 65 | +# # Setup allows N to be solved with SimSolve() given N.groups and target, |
| 66 | +# # golden_section() used to solve outer N.groups component given cost function |
| 67 | +# |
| 68 | +# # Design row conditions must have N and N.groups |
| 69 | +# costfun <- function(condition) |
| 70 | +# with(condition, 5 * N * N.groups + 50 * N.groups) |
| 71 | +# |
| 72 | +# Generate <- function(condition, fixed_objects = NULL){ |
| 73 | +# Attach(condition) |
| 74 | +# group <- rep(letters[1:N.groups], each=N) |
| 75 | +# means <- rnorm(N.groups, sd=.2) |
| 76 | +# DV <- as.vector(sapply(1:N.groups, \(i) rnorm(N, mean=means[i]))) |
| 77 | +# dat <- data.frame(group, DV) |
| 78 | +# dat |
| 79 | +# } |
| 80 | +# |
| 81 | +# Analyse <- function(condition, dat, fixed_objects = NULL) { |
| 82 | +# mod <- aov(DV ~ group, dat) |
| 83 | +# p <- summary(mod)[[1]]$`Pr(>F)`[1L] |
| 84 | +# p |
| 85 | +# } |
| 86 | +# |
| 87 | +# Summarise <- function(condition, results, fixed_objects = NULL) { |
| 88 | +# ret <- EDR(results) |
| 89 | +# ret |
| 90 | +# } |
| 91 | +# |
| 92 | +# \dontrun{ |
| 93 | +# |
| 94 | +# # Solutions given fixed N.groups input |
| 95 | +# Design <- createDesign(N = NA, N.groups=5:25) |
| 96 | +# |
| 97 | +# # test cost function |
| 98 | +# condition <- Design[1,] |
| 99 | +# condition$N <- 200 |
| 100 | +# costfun(condition) |
| 101 | +# |
| 102 | +# # solve for 1-beta = .95 power |
| 103 | +# results <- SimSolve(design=Design, b=.95, interval=c(10, 500), |
| 104 | +# generate=Generate, analyse=Analyse, summarise=Summarise) |
| 105 | +# results |
| 106 | +# |
| 107 | +# # (optional) confirm results accuracy |
| 108 | +# if(FALSE){ |
| 109 | +# cres <- results |
| 110 | +# cres$N <- ceiling(cres$N) |
| 111 | +# confirm <- runSimulation(cres, replications = 10000, |
| 112 | +# generate=Generate, parallel=TRUE, |
| 113 | +# analyse=Analyse, summarise=Summarise) |
| 114 | +# confirm |
| 115 | +# } |
| 116 | +# |
| 117 | +# costs <- numeric(nrow(results)) |
| 118 | +# for(i in 1:length(costs)) |
| 119 | +# costs[i] <- costfun(ceiling(results[i,])) |
| 120 | +# results_costs <- cbind(results, costs) |
| 121 | +# results_costs |
| 122 | +# |
| 123 | +# library(ggplot2) |
| 124 | +# ggplot(results_costs, aes(N, costs, size = N.groups)) + |
| 125 | +# geom_point() |
| 126 | +# |
| 127 | +# |
| 128 | +# # solved via gold_section() with integer option for N.groups |
| 129 | +# Design <- createDesign(N = NA, N.groups = NA) |
| 130 | +# |
| 131 | +# fn <- function(par, condition){ |
| 132 | +# condition$N.groups <- par |
| 133 | +# res <- SimSolve(design=condition, b=.95, interval=c(10, 500), |
| 134 | +# generate=Generate, analyse=Analyse, summarise=Summarise) |
| 135 | +# costfun(res) |
| 136 | +# } |
| 137 | +# |
| 138 | +# res <- gold_section(fn, interval=c(5,25), integer=TRUE, condition=Design) |
| 139 | +# res |
| 140 | +# |
| 141 | +# } |
141 | 142 | gold_section <- function(f, interval, integer = FALSE, |
142 | 143 | tol = .001, maxiter = 100L, ...){ |
143 | 144 | if(integer) interval <- round(interval) |
@@ -166,19 +167,15 @@ gold_section <- function(f, interval, integer = FALSE, |
166 | 167 | xl <- x2; fl <- f2 |
167 | 168 | x2 <- x1; f2 <- f1 |
168 | 169 | x1 <- xl + d.c * (xu - xl) |
169 | | - if(integer){ |
| 170 | + if(integer) |
170 | 171 | x1 <- floor(x1) |
171 | | - if(x1 == x2) break |
172 | | - } |
173 | 172 | f1 <- f(x1, ...) |
174 | 173 | } else { |
175 | 174 | xu <- x1; fu <- f1 |
176 | 175 | x1 <- x2; f1 <- f2 |
177 | 176 | x2 <- xu - d.c * (xu - xl) |
178 | | - if(integer){ |
| 177 | + if(integer) |
179 | 178 | x2 <- ceiling(x2) |
180 | | - if(x1 == x2) break |
181 | | - } |
182 | 179 | f2 <- f(x2, ...) |
183 | 180 | } |
184 | 181 | if(integer) if(abs(xu - xl) <= 1L) break |
|
0 commit comments