MCM(_V)

Description

In the Monte Carlo (MC) method of error propagation, simulated error is added to variables before they are transformed by a user-defined function. This function returns the results of that simulation so that the distribution of the calculated value can be directly observed and tested.

Described here are the function siblings MCM and MCM_V. The former is seedable and nonvolatile while the latter is volatile and uses a slightly better RNG (Excel’s RANDARRAY).

Syntax

MCM(vars, se_vars, fn, [N=1000], [DFs], [skew=0], [kurt], [seed_str]) where

vars is a range containing variable values
se_vars are the standard errors of those variables. There should be as many items in this argument as in vars.
fn is a lambda function that will act on vars and return a scalar value according to the desired transformation of the variables. More details below.
N is the number of simulations to run, defaulting to 1000
DFs is the optional degrees of freedom. If omitted then the normal distribution is used in the simulation (possibly with the Fleishman transformation, if skew or kurtosis is specified). If provided, DFs must be a scalar (in which case it is assumed to apply to all standard errors) or a vector with the same number of items as vars and se_vars. If DFs is provided simulations use the Student’s t-distribution in order to account for uncertainty in the standard errors.
skew and kurt are the desired skewness and excess kurtosis to use in the simulations, obtained by applying the Fleishman method to the transformed normally-distributed RNG. Each can be a single value (used for every variable) or a vector with one value per variable. If skew is given but kurt is omitted, the minimum kurtosis for that skew (rounded up to a whole number; see Fkmin) is used, as with rerror. If DFs is provided these arguments are ignored.
seed_str is the base of a string that is used to generate the random variables, for reproducibility. Since there are multiple calls to the CLS-Stats RNGs in this function, they cannot all use the same seed string so variations are created based on the user-supplied string. The idea is for some reproducibility.

MCM_V has the same syntax as MCM except there is no seed string since it uses Excel’s RANDARRAY RNG, which is not seedable. The result is volatile but a little higher quality and significantly faster than runif.

The lambda function is responsible for extracting the variables from the vars input and then applying the transformation. For example, suppose the desired transformation is

y = a e^{-b}

where a and b are the input variables with some uncertainty $\sigma_a$ and $\sigma_b$ . The vars input would be a row or column vector containing the two input values. The lambda would be defined as (the code block below is formatted for clarity but of course can be all in one line):

LAMBDA(
    vars,			// assuming it is a column vector
    LET(
        a, INDEX(vars,1),
        b, INDEX(vars,2),
        a*EXP(-b)
    )
)

This function is responsible for extracting the individual values from the input range (lines 4-5 above, using Excel’s INDEX function) and then applying the transformation that is being investigated.

If there is only a single variable input then extraction is not necessary.

It might be advisable to create a named lambda function that performs the extraction + transformation since it can be easily tested. Then it would be passed to MCM (without the parentheses) as the function argument, as shown in the example below.

Example

An example of the usage is shown below and in the attached file. In this case, the lambda function was named error_prop and passed by name into the MCM function. Naming the lambda function allowed it to be tested by inputting the values (in cells B2:B3) into the function, which extracted the individual values and applied the $a e^{-b}$ transformation to them. The positive skew of the transformed variables (due to the relatively large relative error on the $b$ variable) is obvious. The MC method will be more accurate than the usual first-order Taylor expansion method commonly used in error propagation. That method also ignores the uncertainty in the standard deviation estimates, but MCM accounts for that by using the t-distribution in the simulation.

MCM for y = a*exp(-b) with a = 100 (SD 5) and b = 1 (SD 0.2): the named lambda returns 36.79 and a QQ plot shows the simulated y is right-skewed

A screenshot of the named lambda function is shown below; the “Refers to” textbox shows the lambda function.

Name Manager showing the error_prop lambda, defined as LAMBDA(vars, LET(a, INDEX(vars,1), b, INDEX(vars,2), a*EXP(-b)))

Of course it could have also been defined in the call, as shown below. Perfectly valid, of course, but this is with only two input variables. With 3+ variables the chances of an error when defining the lambda are greater. The use of the Excel Labs add-in for writing formulas (or defining named lambdas) is recommended in such cases.

=L.MCM(B2:B3,C2:C3,LAMBDA(vars,LET(a,INDEX(vars,1),b,INDEX(vars,2),a*EXP(-b))))

Download: MCM.xlsx

Source Code

Helpers used: FSolve, FSolveBump, FSolveLadder, FSolveMonotonic

The source code for both functions is below. They are very similar except for the fact that the nonvolatile version must handle seeding (using runif) and the volatile function cannot be seeded (and uses RANDARRAY).

MCM = LAMBDA(
    vars, se_vars, fn,
    [N], [DFs], [skew], [kurt], [seed_str],
    LET(
        N, IF(ISOMITTED(N), 1000, N),
        varsRow, TOROW(vars),
        seRow, TOROW(se_vars),
        k, COLUMNS(varsRow),
        useDFs, NOT(ISOMITTED(DFs)),
        seedBase, IF(
            ISOMITTED(seed_str),
            CLS_SystemSeed & ADDRESS(ROW(),COLUMN()) & "MCM",
            seed_str
        ),
        BroadcastRow, LAMBDA(
            x,
            LET(
                v, TOROW(x),
                IF(COLUMNS(v) = 1, MAKEARRAY(1, k, LAMBDA(r,c, INDEX(v,1,1))), v)
            )
        ),
        // length validation
        DFsLenErr, IF(
            NOT(useDFs),
            "",
            IF(
                OR(COLUMNS(TOROW(DFs)) = 1, COLUMNS(TOROW(DFs)) = k),
                "",
                "ERROR: DFs must be a single value or the same length as vars"
            )
        ),
        skewLenErr, IF(
            OR(useDFs, ISOMITTED(skew)),
            "",
            IF(
                OR(COLUMNS(TOROW(skew)) = 1, COLUMNS(TOROW(skew)) = k),
                "",
                "ERROR: skew must be a single value or the same length as vars"
            )
        ),
        kurtLenErr, IF(
            OR(useDFs, ISOMITTED(kurt)),
            "",
            IF(
                OR(COLUMNS(TOROW(kurt)) = 1, COLUMNS(TOROW(kurt)) = k),
                "",
                "ERROR: kurt must be a single value or the same length as vars"
            )
        ),
        lenErr, IF(
            COLUMNS(seRow) <> k,
            "ERROR: vars and se_vars must be the same length",
            IF(DFsLenErr <> "", DFsLenErr, IF(skewLenErr <> "", skewLenErr, kurtLenErr))
        ),
        // skew/kurt prep (only actually evaluated if !useDFs)
        skewRow, IF(ISOMITTED(skew), MAKEARRAY(1, k, LAMBDA(r,c, 0)), BroadcastRow(skew)),
        rawKminRow, Fkmin(skewRow),
        kminRow, IF(rawKminRow < 0, 0, rawKminRow),
        kurtRow, IF(ISOMITTED(kurt), ROUNDUP(kminRow, 0), BroadcastRow(kurt)),
        feasErr, IF(
            useDFs,
            "",
            LET(
                tooSmallRow, kurtRow < kminRow,
                anyTooSmall, OR(tooSmallRow),
                validRow, MAP(
                    SEQUENCE(1,k),
                    LAMBDA(idx, INDEX(FSolveBump(INDEX(skewRow,1,idx), INDEX(kurtRow,1,idx)), 1, 6))
                ),
                anyInvalid, OR(validRow = 0),
                IF(
                    anyTooSmall,
                    "ERROR: kurtosis too small for variable(s) " & TEXTJOIN(",", TRUE, FILTER(SEQUENCE(1,k), tooSmallRow, "")),
                    IF(
                        anyInvalid,
                        "ERROR: no valid distribution found for variable(s) " & TEXTJOIN(",", TRUE, FILTER(SEQUENCE(1,k), validRow = 0, "")),
                        ""
                    )
                )
            )
        ),
        // DFs prep (only actually evaluated if useDFs)
        dfRow, IF(useDFs, BroadcastRow(DFs), MAKEARRAY(1, k, LAMBDA(r,c, 0))),
        allErr, IF(lenErr <> "", lenErr, feasErr),
        IF(
            allErr <> "",
            allErr,
            LET(
                firstCol, IF(
                    useDFs,
                    T.INV(runif(N,1,,,,"v_" & 1 & seedBase), INDEX(dfRow,1,1)),
                    rerror(N, 1, INDEX(skewRow,1,1), INDEX(kurtRow,1,1), "v_" & 1 & seedBase)
                ),
                noiseMatrix, IF(
                    k = 1,
                    firstCol,
                    REDUCE(
                        firstCol,
                        SEQUENCE(1, k-1),
                        LAMBDA(acc, idx,
                            LET(
                                j, idx + 1,
                                colJ, IF(
                                    useDFs,
                                    T.INV(runif(N,1,,,,"v_" & j & seedBase), INDEX(dfRow,1,j)),
                                    rerror(N, 1, INDEX(skewRow,1,j), INDEX(kurtRow,1,j), "v_" & j & seedBase)
                                ),
                                HSTACK(acc, colJ)
                            )
                        )
                    )
                ),
                simMatrix, varsRow + seRow * noiseMatrix,
                BYROW(simMatrix, fn)
            )
        )
    )
);

MCM_V = LAMBDA(
    vars, se_vars, fn,
    [N], [DFs], [skew], [kurt],
    LET(
        N, IF(ISOMITTED(N), 1000, N),
        varsRow, TOROW(vars),
        seRow, TOROW(se_vars),
        k, COLUMNS(varsRow),
        useDFs, NOT(ISOMITTED(DFs)),
        BroadcastRow, LAMBDA(
            x,
            LET(
                v, TOROW(x),
                IF(COLUMNS(v) = 1, MAKEARRAY(1, k, LAMBDA(r,c, INDEX(v,1,1))), v)
            )
        ),
        DFsLenErr, IF(
            NOT(useDFs),
            "",
            IF(
                OR(COLUMNS(TOROW(DFs)) = 1, COLUMNS(TOROW(DFs)) = k),
                "",
                "ERROR: DFs must be a single value or the same length as vars"
            )
        ),
        skewLenErr, IF(
            OR(useDFs, ISOMITTED(skew)),
            "",
            IF(
                OR(COLUMNS(TOROW(skew)) = 1, COLUMNS(TOROW(skew)) = k),
                "",
                "ERROR: skew must be a single value or the same length as vars"
            )
        ),
        kurtLenErr, IF(
            OR(useDFs, ISOMITTED(kurt)),
            "",
            IF(
                OR(COLUMNS(TOROW(kurt)) = 1, COLUMNS(TOROW(kurt)) = k),
                "",
                "ERROR: kurt must be a single value or the same length as vars"
            )
        ),
        lenErr, IF(
            COLUMNS(seRow) <> k,
            "ERROR: vars and se_vars must be the same length",
            IF(DFsLenErr <> "", DFsLenErr, IF(skewLenErr <> "", skewLenErr, kurtLenErr))
        ),
        skewRow, IF(ISOMITTED(skew), MAKEARRAY(1, k, LAMBDA(r,c, 0)), BroadcastRow(skew)),
        rawKminRow, Fkmin(skewRow),
        kminRow, IF(rawKminRow < 0, 0, rawKminRow),
        kurtRow, IF(ISOMITTED(kurt), ROUNDUP(kminRow, 0), BroadcastRow(kurt)),
        feasErr, IF(
            useDFs,
            "",
            LET(
                tooSmallRow, kurtRow < kminRow,
                anyTooSmall, OR(tooSmallRow),
                validRow, MAP(
                    SEQUENCE(1,k),
                    LAMBDA(idx, INDEX(FSolveBump(INDEX(skewRow,1,idx), INDEX(kurtRow,1,idx)), 1, 6))
                ),
                anyInvalid, OR(validRow = 0),
                IF(
                    anyTooSmall,
                    "ERROR: kurtosis too small for variable(s) " & TEXTJOIN(",", TRUE, FILTER(SEQUENCE(1,k), tooSmallRow, "")),
                    IF(
                        anyInvalid,
                        "ERROR: no valid distribution found for variable(s) " & TEXTJOIN(",", TRUE, FILTER(SEQUENCE(1,k), validRow = 0, "")),
                        ""
                    )
                )
            )
        ),
        dfRow, IF(useDFs, BroadcastRow(DFs), MAKEARRAY(1, k, LAMBDA(r,c, 0))),
        allErr, IF(lenErr <> "", lenErr, feasErr),
        IF(
            allErr <> "",
            allErr,
            LET(
                firstCol, IF(
                    useDFs,
                    T.INV(RANDARRAY(N,1), INDEX(dfRow,1,1)),
                    rerrorV(N, 1, INDEX(skewRow,1,1), INDEX(kurtRow,1,1))
                ),
                noiseMatrix, IF(
                    k = 1,
                    firstCol,
                    REDUCE(
                        firstCol,
                        SEQUENCE(1, k-1),
                        LAMBDA(acc, idx,
                            LET(
                                j, idx + 1,
                                colJ, IF(
                                    useDFs,
                                    T.INV(RANDARRAY(N,1), INDEX(dfRow,1,j)),
                                    rerrorV(N, 1, INDEX(skewRow,1,j), INDEX(kurtRow,1,j))
                                ),
                                HSTACK(acc, colJ)
                            )
                        )
                    )
                ),
                simMatrix, varsRow + seRow * noiseMatrix,
                BYROW(simMatrix, fn)
            )
        )
    )
);

Fkmin MCM_CI(_V)