monteCarloSimulation {stepR}

Monte Carlo simulation

Description

Performs Monte-Carlo simulations of the multiscale vector of statistics, (3.9) in the vignette, and of the penalised multiscale statistic, (3.6) in the vignette, when no signal is present, see also section 3.2.3 in the vignette.

Usage

monteCarloSimulation(n, r = 1e4L, family = NULL, intervalSystem = NULL,
                     lengths = NULL, penalty = NULL,
                     output = c("vector", "maximum"), seed = n,
                     rand.gen = NULL, messages = NULL, ...)

Arguments

n

a positive integer giving the number of observations for which the Monte-Carlo simulation will be performed

r

a positive integer giving the number of repititions

family

a string specifying the assumed parametric family, for more details see parametricFamily, currently "gauss", "hsmuce" and "mDependentPS" are supported. By default (NULL) "gauss" is assumed

intervalSystem

a string giving the used interval system, either "all" for all intervals, "dyaLen" for all intervals of dyadic length or "dyaPar" for the dyadic partition, for more details see intervalSystem. By default (NULL) the default interval system of the specified parametric family will be used, which one this will be is described in parametricFamily

lengths

an integer vector giving the set of lengths, i.e. only intervals of these lengths will be considered. Only required for output == "maximum", otherwise ignored with a warning. Note that not all lengths are possible for all interval systems and for all parametric families, see intervalSystem and parametricFamily, respectively, to see which ones are allowed. By default (NULL) all lengths that are possible for the specified intervalSystem and for the specified parametric family will be used

penalty

a string specifying how the statistics will be penalised, either "sqrt", "log" or "none", see penalty and section 3.2 in the vignette for more details. Only required for output == "maximum", otherwise ignored with a warning. By default (NULL) the default penalty of the specified parametric family will be used, which one this will be is described in parametricFamily

output

a string specifying the output, see Value

seed

will be passed to set.seed to set a seed, set.seed will not be called if this argument is set to "no", i.e. a single value, interpreted as an integer, NULL or "no"

rand.gen

by default (NULL) this argument will be replaced by the default function to generate random observations of the given family. Note that a Monte-Carlo simulation can only be saved if rand.gen == NULL. Alternatively, an own function expecting a single argument named data and returning a numeric vector of length n, this is given by data$n. Will be called with rand.gen(data = data), with data a list containing the named entries n, the expected number of data points, and parameters of the parametric family, e.g. sd for family == "gauss" or covariances for family == "mDependentPS"

messages

a positive integer or NULL, in each messages iteration a message will be printed in order to show the progress of the simulation, if NULL no message will be given

...

further parameters of the parametric family. Depending on the argument family some might be required, but others might be optional, please see parametricFamily for more details

Value

If output == "vector" an object of class "MCSimulationVector", i.e. a d_n times r matrix containing r independent samples of the multiscale vector of statistics, with d_n the number of scales, i.e. the number of possible lengths for the given interval system and given parametric family. If output == "maximum" an object of class "MCSimulationMaximum", i.e. a vector of length r containing r independent samples of the penalised multiscale statistic. For both, additionally, the following attributes are set:

• "keyList": A list specifying for which number of observations n, which parametric family with which parameters by a SHA-1 hash, which interval system and in case of "MCSimulationMaximum", additionally, for which lengths and which penalisation the simulation was performed.

• "key": A key used internally for identification when the object will be saved and loaded.

• "n": The number of observations n for which the simulation was performed.

• "lengths": The lengths for which the simulation was performed.

• "save": A logical which is TRUE if the object can be saved which is the case for rand.gen == NULL and FALSE otherwise.

References

Frick, K., Munk, A., Sieling, H. (2014) Multiscale change-point inference. With discussion and rejoinder by the authors. Journal of the Royal Statistical Society, Series B 76(3), 495–580.

Pein, F., Sieling, H., Munk, A. (2017) Heterogeneous change point inference. Journal of the Royal Statistical Society, Series B, 79(4), 1207–1227.

See Also

critVal, computeStat, penalty, parametricFamily, intervalSystem

Examples

# monteCarloSimulation will be called in critVal, can be called explicitly
# object of class MCSimulationVector
stat <- monteCarloSimulation(n = 100L)

identical(critVal(n = 100L, alpha = 0.5, stat = stat),
          critVal(n = 100L, alpha = 0.5,
                  options = list(load = list(), simulation = "matrix")))

# object of class MCSimulationMaximum
stat <- monteCarloSimulation(n = 100L, output = "maximum")
identical(critVal(n = 100L, alpha = 0.5, stat = stat),
          critVal(n = 100L, alpha = 0.5,
                  options = list(load = list(), simulation = "vector")))

# different interval system, lengths and penalty
monteCarloSimulation(n = 100L, output = "maximum", intervalSystem = "dyaLen",
              lengths = c(1L, 2L, 4L, 8L), penalty = "log")

# with a different number of iterations, different seed,
# reported progress and user written rand.gen function
stat <- monteCarloSimulation(n = 100L, r = 1e3, seed = 1, messages = 100,
                             rand.gen = function(data) {rnorm(100)})

# the optional argument sd of parametric family "gauss" will be replaced by 1
identical(monteCarloSimulation(n = 100L, r = 1e3, sd = 5),
          monteCarloSimulation(n = 100L, r = 1e3, sd = 1))

# simulation for family "hsmuce"
monteCarloSimulation(n = 100L, family = "hsmuce")

# simulation for family "mDependentGauss"
# covariances must be given (can also be given by correlations or filter)
stat <- monteCarloSimulation(n = 100L, family = "mDependentPS",
                             covariances = c(1, 0.5, 0.3))

# variance will be standardized to 1
# output might be on some systems even identical
all.equal(monteCarloSimulation(n = 100L, family = "mDependentPS",
                               covariances = c(2, 1, 0.6)), stat)