# Monte Carlo simulation

I have an equation to calculate `z` from multiplying several rasters and numbers as seen below:

``````z <- a*b*c*d*e
``````

The raster layers and values are as below:

``````library(terra)
set.seed(234)
r <- rast(nrows=10, ncols=10)

a <- setValues(r,runif(n = 100, min = 4, max=16))
b <- setValues(r,runif(n = 100, min = 0, max=100))
c <- setValues(r,runif(n = 100, min = 68, max=270))
d <- 0.8
e <- 7.9
``````

I need to calculate uncertainty based on varying `c` and `e` by varying both variables as below:

`+/- 10` for `c`

`+/-` 4 for `e`

values need to change independently. In the end, I need to plot an uncertainty map.

• Have you tried? Write a function that computes `z` given `c` and `e` and then loop over random realisations of those, keeping stats. Although since `z` is pretty simply defined you can probably work out E(z) and Var(z) given some distribution of `c` and `e` with the tiniest bit of statistical theory and avoid the Monte-Carlo. Commented Sep 12, 2022 at 11:29
• @Spacedman could you share a sample of how the process would look like? Commented Sep 12, 2022 at 11:31
• What distribution for `c` and `e` do you have? Are they `Uniform` between those limits? Or Normally distributed with mean 0 and given SD? Are they independent or correlated? Commented Sep 12, 2022 at 11:34
• They both have normal distribution Commented Sep 12, 2022 at 11:36
• The range of a Normal distribution is `-Inf, +Inf`, so you should edit your question and say what the variance or standard deviation of the quantities are. Commented Sep 12, 2022 at 11:39

The below may give you some ideas. I am using normal distributions with your +/- as standard deviations. I assumed that you want the values of all cells of `c` to change in the same way in each iteration, but you may want them to change independently instead.

``````library(terra)
set.seed(234)
r <- rast(nrows=10, ncols=10)

a <- setValues(r,runif(n = 100, min = 4, max=16))
b <- setValues(r,runif(n = 100, min = 0, max=100))
c <- setValues(r,runif(n = 100, min = 68, max=270))
d <- 0.8
e <- 7.9
``````

Solution

``````f <- function(n) {
ee <- rnorm(n, e, 4)
dc <- rnorm(n, 0, 10)
(a * b * d) * (c + dc) * ee
}

s <- f(100)
q <- quantile(s)
q
#class       : SpatRaster
#dimensions  : 10, 10, 5  (nrow, ncol, nlyr)
#resolution  : 36, 18  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84
#source      : memory
#names       :         q0,     q0.25,       q0.5,      q0.75,         q1
#min values  : -636983.52,   22728.56,   32998.79,   41095.48,   79844.2
#max values  :  -11830.85, 1319820.92, 1823798.55, 2304815.29, 4404378.6

plot(q)
``````