/
distribution.R
437 lines (385 loc) · 13.9 KB
/
distribution.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
#_______________________________________________________________________________
#---- distribution class ----
#_______________________________________________________________________________
validateDistribution <- function(object) {
return(TRUE)
}
#'
#' Distribution class. See this class as an interface.
#'
#' @export
setClass(
"distribution",
representation(
),
contains="pmx_element",
validity=validateDistribution
)
#_______________________________________________________________________________
#---- undefined_distribution class ----
#_______________________________________________________________________________
#'
#' Undefined distribution class. This type of object is automatically created
#' in method toExplicitDistribution() when the user does not provide a concrete
#' distribution. This is because S4 objects do not accept NULL values.
#'
#' @export
setClass(
"undefined_distribution",
representation(
),
contains="distribution"
)
#_______________________________________________________________________________
#---- sampled_distribution class ----
#_______________________________________________________________________________
validateSampledDistribution <- function(object) {
return(expectZeroOrMore(object, "sampled_values"))
}
setClass(
"sampled_distribution",
representation(
sampled_values = "numeric"
),
contains="distribution",
prototype=prototype(sampled_values=numeric(0)),
validity=validateSampledDistribution
)
#_______________________________________________________________________________
#---- constant_distribution class ----
#_______________________________________________________________________________
validateConstantDistribution <- function(object) {
return(expectOneForAll(object, c("value")))
}
#'
#' Constant distribution class.
#'
#' @slot value covariate value, single numeric value
#' @export
setClass(
"constant_distribution",
representation(
value = "numeric"
),
contains="sampled_distribution",
validity=validateConstantDistribution
)
#'
#' Create a constant distribution. Its value will be constant across all generated samples.
#'
#' @param value covariate value, single numeric value
#' @return a constant distribution (same value for all samples)
#' @export
ConstantDistribution <- function(value) {
return(new("constant_distribution", value=as.numeric(value)))
}
#_______________________________________________________________________________
#---- fixed_distribution class ----
#_______________________________________________________________________________
validateFixedDistribution <- function(object) {
return(expectOneOrMore(object, c("values")))
}
#'
#' Fixed distribution class.
#'
#' @slot values covariate values, numeric vector (1 value per sample)
#' @export
setClass(
"fixed_distribution",
representation(
values = "numeric"
),
contains="sampled_distribution",
validity=validateFixedDistribution
)
#'
#' Create a fixed distribution.
#' Each sample will be assigned a fixed value coming from vector 'values'.
#'
#' @param values covariate values, numeric vector (1 value per sample)
#' @return a fixed distribution (1 value per sample)
#' @export
FixedDistribution <- function(values) {
return(new("fixed_distribution", values=values))
}
#_______________________________________________________________________________
#---- function_distribution class ----
#_______________________________________________________________________________
validateFunctionDistribution <- function(object) {
check1 <- expectOne(object, "fun")
check2 <- expectZeroOrMore(object, "args")
return(c(check1, check2))
}
#'
#' Function distribution class.
#'
#' @slot fun function name, character (e.g. 'rnorm')
#' @slot args list of arguments (e.g list(mean=70, sd=10))
#' @export
setClass(
"function_distribution",
representation(
fun = "character",
args = "list"
),
contains="sampled_distribution",
prototype=prototype(args=list()),
validity=validateFunctionDistribution
)
#'
#' Create a function distribution. During distribution sampling, the provided function
#' will be responsible for generating values for each sample. If first argument
#' of this function is not the size (n), please tell which argument corresponds
#' to the size 'n' (e.g. list(size="n")).
#'
#' @param fun function name, character (e.g. 'rnorm')
#' @param args list of arguments (e.g list(mean=70, sd=10))
#' @return a function distribution
#' @export
FunctionDistribution <- function(fun, args) {
return(new("function_distribution", fun=fun, args=args))
}
#'
#' Create an uniform distribution.
#'
#' @param min min value
#' @param max max value
#' @return an uniform distribution
#' @export
UniformDistribution <- function(min, max) {
expectSingleNumericValue(min, "min")
expectSingleNumericValue(max, "max")
return(new("function_distribution", fun="runif", args=list(min=as.numeric(min), max=as.numeric(max))))
}
#'
#' Create a normal distribution.
#'
#' @param mean mean value of distribution
#' @param sd standard deviation of distribution
#' @return a normal distribution
#' @export
NormalDistribution <- function(mean, sd) {
expectSingleNumericValue(mean, "mean")
expectSingleNumericValue(sd, "sd")
return(new("function_distribution", fun="rnorm", args=list(mean=as.numeric(mean), sd=as.numeric(sd))))
}
#'
#' Create a log normal distribution.
#'
#' @param meanlog mean value of distribution in log domain
#' @param sdlog standard deviation of distribution in log domain
#' @return a log normal distribution
#' @export
LogNormalDistribution <- function(meanlog, sdlog) {
expectSingleNumericValue(meanlog, "meanlog")
expectSingleNumericValue(sdlog, "sdlog")
return(new("function_distribution", fun="rlnorm", args=list(meanlog=as.numeric(meanlog), sdlog=as.numeric(sdlog))))
}
#'
#' Discrete distribution.
#'
#' @param x vector of one or more integers from which to choose
#' @param prob a vector of probability weights for obtaining the elements of the vector being sampled
#' @param replace should sampling be with replacement, default is TRUE
#' @return a discrete distribution
#' @export
DiscreteDistribution <- function(x, prob, replace=TRUE) {
assertthat::assert_that(is.numeric(x), msg="x must be numeric")
assertthat::assert_that(is.numeric(prob), msg="prob must be numeric")
assertthat::assert_that(length(x)==length(prob), msg="x and prob must have the same length")
return(new("function_distribution", fun="base::sample", args=list(size="n", x=as.numeric(x), prob=as.numeric(prob), replace=as.logical(replace))))
}
#'
#' Binomial distribution.
#'
#' @param trials number of Bernoulli trials per observation (=subject), integer
#' @param prob probability of success for each trial
#' @return a binomial distribution
#' @export
BinomialDistribution <- function(trials, prob) {
expectSingleIntegerValue(trials, name="trials")
expectSingleNumericValue(prob, name="prob")
return(new("function_distribution", fun="rbinom", args=list(size=as.integer(trials), prob=as.numeric(prob))))
}
#'
#' Retrieve the parameter value (standardized) for the specified parameter name.
#'
#' @param model model
#' @param paramName parameter name
#' @param default default value if not found
#' @param mandatory must be in model or not
#' @return the standardized parameter value or the given default value if not found
#' @importFrom assertthat assert_that
#' @export
retrieveParameterValue <- function(model, paramName, default=NULL, mandatory=FALSE) {
assertthat::assert_that(is.character(paramName) && length(paramName)==1,
msg="paramName must be a single character value")
parameter <- model@parameters %>% getByName(paramName)
if (parameter %>% length() == 0) {
if (is.null(default) && mandatory) {
stop(paste0(paramName, " can't be found in model"))
}
return(default)
} else {
# If parameter is OMEGA, it needs to be standardised before taking its value
# This way, value is always a variance (or covariance)
parameter <- parameter %>% standardise()
return(parameter@value)
}
}
#'
#' Create a parameter distribution. The resulting distribution is a
#' log-normal distribution, with meanlog=log(THETA) and sdlog=sqrt(OMEGA).
#'
#' @param model model
#' @param theta corresponding THETA name, character
#' @param omega corresponding OMEGA name, character, NULL if not defined
#' @return a parameter distribution
#' @export
ParameterDistribution <- function(model, theta, omega=NULL) {
thetaValue <- retrieveParameterValue(model, paramName=paste0("THETA_", theta), mandatory=TRUE)
if (is.null(omega)) {
omegaValue <- 0
} else {
omegaValue <- retrieveParameterValue(model, paramName=paste0("OMEGA_", omega), mandatory=TRUE)
}
fun <- FunctionDistribution(fun="rlnorm", args=list(meanlog=log(thetaValue), sdlog=sqrt(omegaValue)))
return(fun)
}
#'
#' Create an ETA distribution. The resulting distribution is a
#' normal distribution, with mean=0 and sd=sqrt(OMEGA).
#'
#' @param model model
#' @param omega corresponding THETA name, character
#' @return an ETA distribution
#' @export
EtaDistribution <- function(model, omega) {
omegaValue <- retrieveParameterValue(model, paramName=paste0("OMEGA_", omega), mandatory=TRUE)
return(NormalDistribution(mean=0, sd=sqrt(omegaValue)))
}
#_______________________________________________________________________________
#---- bootstrap_distribution class ----
#_______________________________________________________________________________
validateBootstrapDistribution <- function(object) {
check1 <- expectOneOrMore(object, c("data"))
check2 <- expectOneForAll(object, c("replacement", "random"))
return(c(check1, check2))
}
#'
#' Bootstrap distribution class.
#'
#' @slot data values to draw, numeric vector
#' @slot replacement values can be reused or not, logical
#' @slot random values are drawn randomly, logical
#' @export
setClass(
"bootstrap_distribution",
representation(
data = "numeric",
replacement = "logical",
random = "logical"
),
contains="sampled_distribution",
prototype=prototype(replacement=FALSE, random=FALSE),
validity=validateBootstrapDistribution
)
#'
#' Create a bootstrap distribution. During function sampling, CAMPSIS will generate
#' values depending on the given data and arguments.
#'
#' @param data values to draw, numeric vector
#' @param replacement values can be reused or not, logical
#' @param random values are drawn randomly, logical
#' @return a bootstrap distribution
#' @export
BootstrapDistribution <- function(data, replacement=FALSE, random=FALSE) {
return(new("bootstrap_distribution", data=data, replacement=replacement, random=random))
}
#_______________________________________________________________________________
#---- sample ----
#_______________________________________________________________________________
sample_delegate <- function(fun, ...) {
return(fun(...))
}
sample_delegate_n <- function(fun, n, ...) {
return(fun(n, ...))
}
#' @rdname sample
setMethod("sample", signature = c("constant_distribution", "integer"), definition = function(object, n) {
object@sampled_values <- rep(object@value, n)
return(object)
})
#' @rdname sample
setMethod("sample", signature = c("fixed_distribution", "integer"), definition = function(object, n) {
object@sampled_values <- object@values
if (length(object@values) != n) {
stop(paste0("A fixed distribution should have exactly ", n, " values, not ", length(object@values)))
}
return(object)
})
#' @rdname sample
setMethod("sample", signature = c("function_distribution", "integer"), definition = function(object, n) {
fun <- eval(parse(text=object@fun))
args <- object@args
SAMPLING_SIZE <- n
# If "n" is provided as value in args, this means the first arg of the function
# is not related to the size
hasNCharAsValue <- "n" %in% as.character(args)
if (length(args) == 0) {
args_str <- "SAMPLING_SIZE"
} else {
# Example: return 'mean=70, sd=2'
args_str <- purrr::accumulate2(names(args), args, .f=function(.x, .y, .z){
comma <- if (.x=="") {""} else {", "}
if (.z[[1]]=="n") {
.z <- "SAMPLING_SIZE"
}
if (length(.z) > 1) {
value <- paste0("c(", paste0(.z, collapse=","), ")")
} else {
value <- .z
}
paste0(.x, comma, .y , "=", value)
}, .init="")
if (hasNCharAsValue) {
# Size is not the first argument
args_str <- args_str[length(args_str)]
} else {
# Size is the first argument
args_str <- paste0("SAMPLING_SIZE, ", args_str[length(args_str)])
}
}
# Eval string expression
if (hasNCharAsValue) {
text <- paste0("sample_delegate(fun, ", args_str, ")")
} else {
text <- paste0("sample_delegate_n(fun, ", args_str, ")")
}
values <- eval(parse(text=text))
# Check result
if (length(values) != n) {
stop(paste("Calling function", object@fun, "with arguments", args_str, "does not provided", n, "values but", length(values)))
}
# Assign values
object@sampled_values <- values
return(object)
})
#' @rdname sample
setMethod("sample", signature = c("bootstrap_distribution", "integer"), definition = function(object, n) {
data <- object@data
nData <- length(data)
replacement <- object@replacement
random <- object@random
if (n > nData && !replacement) {
stop(paste("Provided bootstrap function does not have enough data to generate", n, "samples (please set 'replacement' to TRUE)"))
}
if (random) {
values <- data[sample.int(n=nData, size=n, replace=replacement)]
} else {
values <- rep(data, length.out=n)
}
# Assign values
object@sampled_values <- values
return(object)
})