Monte Carlo simulation

Question

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.

Answer 1

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.

Your example data

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)