MC.analysis: Analysis of the Monte Carlo simulation

Description Usage Arguments Value Author(s) See Also Examples

View source: R/MC.analysis.R

Description

Function for running the analysis of the Monte Carlo simulation.

Usage

1
2
  MC.analysis(x, delta, qUpper, p1.det, sim.det, event.ini, event.end, 
  ntick, summ.data = NULL)

Arguments

x

A list .

delta

A numeric value that specifies the level of aggregation required in minutes.

qUpper

A character string that defines the upper percentile to plot the confidence band of results, several options are possible "q999" the 99.9th percentile, "q995" the 99.5th percentile, "q99" the 99th percentile, "q95" the 95th percentile, "q50" the 50th percentile. The lower boundary of the confidence band (showed in gray in the output plots) is the 5th percentile in all cases.

p1.det

A data.frame that contains the time series of the main driving force of the system to be simulated deterministically, e.g. precipitation. This data.frame should have only two columns: the first one, Time [y-m-d h:m:s]; the second one, a numeric value equal to the magnitude of the environmental variable.

sim.det

A list that contains the results of the deterministic simulation, here the output of EmiStatR given p1.det. See the method EmiStatR from the homonym package for details.

event.ini

A time-date string in POSIXct format that defines the initial time for event analysis.

event.end

A time-date string in POSIXct format that defines the final time for event analysis.

ntick

A numeric value to specify the number of ticks in the x-axis for the event time-window plots.

summ.data

A list by default NULL. If provided, the list should contain an output of the MC.analysis function, and the analysis is done again without the calculation of some of the internal variables, therefore the analysis is faster.

Value

A list of length 2:

summ

A list that contains the summary statistics of the Monte Carlo simulation per output variable. Each output variable is summarised by calculating the mean "Mean", standard deviation "sd", variance "Variance", 5th, 25th, 50th, 75th, 95th, 99.5th, 99.9th percentiles "q05", "q25", "q50", "q75", "q95", "q995", "q999", the max "Max", the sum "Sum", time "time", and the deterministic precipitation "p1", all variables as time series.

variance

A data.frame that contains the summary statistics of the variance of the Monte Carlo simulation per output variable.

Author(s)

J.A. Torres-Matallana

See Also

See also setup-class, MC.setup-methods, MC.sim-methods.

Examples

  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
## the Monte Carlo simulation: MC.sim
library(EmiStatR)

# library(xts)
# data(Esch_Sure2010)
# P <- IsReg(Esch_Sure2010, format = "%Y-%m-%d %H:%M:%S", tz = "CET")
# P1 <- P[[2]]
# P1 <- P1["2010-08",][1:55]
# P1 <- cbind.data.frame(time=index(P1), P1 = coredata(P1))

data(P1)
P1 <- P1[165:(110*2),]
plot(P1[,2], typ="l")

library(stUPscales)

setting_EmiStatR <-  setup(id       = "MC_sim1",
                           nsim     = 3, # # use a larger number to have 
                                           # a proper confidence band of simulatios
                           seed     = 123, 
                           mcCores  = 1, 
                           ts.input = P1,  
                           rng      = rng <- list(
                             qs   = 150,    # [l/PE/d]
                             CODs = c(pdf = "nor", mu = 4.378, sigma = 0.751),    # log[g/PE/d]
                             NH4s = c(pdf = "nor", mu = 1.473, sigma = 0.410),    # log[g/PE/d]
                             qf   = 0.04,     # [l/s/ha]
                             CODf = 0,              # [g/PE/d]
                             NH4f = 0,              # [g/PE/d]
                             CODr = c(pdf = "nor", mu = 3.60, sigma = 1.45),     # 71 log[mg/l]
                             NH4r = 1,              # [mg/l]
                             nameCSO = "E1",        # [-]
                             id      = 1,           # [-]
                             ns      = "FBH Goesdorf",  # [-]
                             nm      = "Goesdorf",  # [-]
                             nc      = "Obersauer", # [-]
                             numc    = 1,          # [-]
                             use     = "R/I",  # [-]
                             Atotal  = 36,              # [ha]
                             Aimp    = c(pdf = "uni", min = 4.5, max = 25),       # [ha]
                             Cimp    = c(pdf = "uni", min = 0.25, max = 0.95),  # [-]
                             Cper    = c(pdf = "uni", min = 0.05, max = 0.60),  # [-]
                             tfS     = 1,               # [time steps]
                             pe      = 650,             # [PE]
                             Qd      = 5,               # [l/s]
                             Dd      = 0.150,           # [m]
                             Cd      = 0.18,            # [-]
                             V       = 190,             #  [m3]
                             lev.ini = 0.10,            # [m]
                             lev2vol = list(lev = c(.06, 1.10, 1.30, 3.30),   # [m]
                                            vol = c(0, 31, 45, 190))          # [m3]
                           ),
                           ar.model  = ar.model <- list(
                             CODs    = 0.5,         
                             NH4s    = 0.5,
                             CODr    = 0.7),
                           var.model = var.model <- list(
                             inp     = c("", ""), # c("CODs", "NH4s"), # c("", ""),
                             w       = c(0.04778205, 0.02079010),
                             A       = matrix(c(9.916452e-01, -8.755558e-05, 
                                                -0.003189094, 0.994553910), nrow=2, ncol=2),
                             C       = matrix(c(0.009126591, 0.002237936, 
                                                0.002237936, 0.001850941), nrow=2, ncol=2)))

MC_setup <- MC.setup(setting_EmiStatR)

sims <- MC.sim(x = MC_setup, EmiStatR.cores = 0)


## Monte Carlo simulation analysis: MC.analysis

# Deterministic simulation
# Definition of structure 1, E1:

E1 <- list(id = 1, ns = "FBH Goesdorf", nm = "Goesdorf", nc = "Obersauer", numc = 1,
           use = "R/I", Atotal = 36, Aimp = 25.2, Cimp = 0.80, Cper = 0.30,
           tfS = 0, pe = 650, Qd = 5, 
           Dd = 0.150, Cd = 0.18, V = 190, lev.ini = 0.10,
           lev2vol = list(lev = c(.06, 1.10, 1.30, 3.30),   
                          vol = c(0, 31, 45, 190))
           )

# Defining deterministic input:
library(EmiStatR)
# data(P1)

input.det <- input(spatial = 0, zero = 1e-5, 
                    folder = system.file("shiny", package = "EmiStatR"),
                    cores = 0,
                    ww = list(qs = 150, CODs = 104, NH4s = 4.7),
                    inf = list(qf= 0.04, CODf = 0, NH4f = 0),
                    rw = list(CODr = 71, NH4r = 1, stat = "Dahl"),
                    P1 = P1, st = list(E1=E1), export = 0)

# Invoking `EmiStatR` with the deterministic input: 
sim.det   <- EmiStatR(input.det)

# further arguments
delta <- 10 # minutes 
qUpper <- "q999" 

event.ini <- as.POSIXct("2016-01-02 03:20:00")
event.end <- as.POSIXct("2016-01-02 12:30:00")

# uncomment to run:
# new_analysis <- MC.analysis(x = sims, delta = delta, qUpper = qUpper,  p1.det = P1, 
#                             sim.det = sim.det, event.ini = event.ini, event.end = event.end, 
#                             ntick = 5, summ.data = NULL)

stUPscales documentation built on May 6, 2019, 1:08 a.m.