How to find variance of multivariable expression I have an expression of this form
$W_i = a * X_i * (1 - Y_i * (1 - Z_i))$,
where $a$ is a constant, $X_i \sim Bernoulli({p_x}_i)$,
$Y_i \sim Bernoulli({p_y}_i)$, and
$Z_i \sim N(\mu_i, \sigma)$
Background:
So, I have three models for X, Y and Z, and these are models for probability of survival, probability of death, and a model for some distance measure, respectively. So X and Y are logistic regression models, and Z is a normal regression model. I have a dataset, and I predict the expected values for the three models, and then I calculate W. But I know that this expected W has a high degree of uncertainty attached to it, because the underlying drivers are highly uncertain. I wish to get a distribution of W for all datapoints, so I can measure the uncertainty.
Attempt:
I have 100k rows of data, so $i = {1,2,3,...,100000}$, where I have predicted values in each row for ${p_x}_i$, ${p_y}_i$ and $\mu_i$, and $\sigma$ is constant/equal for all rows. I want to find the variance of the expression for W above for each row, and I've tried Monte Carlo sampling, and then bootstrapping to find the variance, but it is just very slow.
Is there a faster way to do this, or do I need to calculate $Var(W_i)$ "by hand"?
Code:
# data
set.seed(1234)
n <- 100000
p_x <- rnorm(n, 0.3, 0.02)
p_y <- rnorm(n, 0.7, 0.02)
mu_z <- rnorm(n, 0.5, 0.05)
sigma = 0.1
a = 10

# monte carlo simulations
## duplicate each i/observation n_mc times using rep(), and sample
n_mc = 100
X <- rbinom(n * n_mc, 1, rep(p_x, each = n_mc))
Y <- rbinom(n * n_mc, 1, rep(p_y, each = n_mc))
Z <- rnorm(n * n_mc, rep(mu_z, each = n_mc), sd = rep(sigma, n * n_mc))
W <- X * (1 - (Y * (1 - Z))) * a
W <- matrix(W, nrow=n, ncol=n_mc, byrow=TRUE)

# bootstrap
## since each row in W is a sample, we bootstrap this sample to obtain variance/or confidence intervals
library(boot)
getMean <- function(x, ind){
  mean(x[ind], na.rm = TRUE)
}
out <- do.call(rbind, apply(W, 1, FUN = function(x) 
      {
    x1 <- boot.ci(boot(x, statistic = getMean, R = 1000), conf = 0.8, type="perc")
    data.frame(mean = x1$t0, lowerCI = x1$percent[4], upperCI = x1$percent[5])
      }))


This is a partial answer and partially a request for clarification.
If I understand correctly,

*

*$a$ is a constant.


*$X_i$ is Bernoulli with $\Pr(X_i=1)=p_i.$


*$Y_i$ is Bernoulli with $\Pr(Y_i=1)=q_i.$  (I changed the names of the parameters slightly to avoid cascading subscripts.)


*$Z_i$ has a Normal$(\mu_i,\sigma)$ distribution.


*All the variables are independent.  (If not, we don't have enough information to solve this problem.)
To evaluate the distribution of $W_i=a X_i(1 - Y_i(1-Z_i)),$ ignore $a$ for the moment because all it does is scale everything and set $U_i = X_i(1 - Y_i(1-Z_i)).$  We will apply definitions and basic properties, beginning by noting that $(X_i,Y_i)$ is certain to be one of just four values, as presented in this table.
$$\begin{array}{rrlll}
& \Pr & X_i & Y_i & U_i\\ \hline
& (1-p_i)(1-q_i) & 0 & 0 & 0\\
& (1-p_i)q_i & 0 & 1 & 0\\
& p_i(1-q_i) & 1 & 0 & 1\\
& p_iq_i & 1 & 1 & Z_i\\ \hline
\end{array}$$
This evidently is a mixture of an atom at $0$ with probability $(1-p_i)(1-q_i)+(1-p_i)q_i = q_i,$ an atom at $1$ with probability $p_i(1-q_i),$ and a Normal variable $V_i.$  The calculation of its moments is immediate, so that in particular
$$E[U_i] = q_i(0) + p_i(1-q_i)(1) + p_iq_iE[Z_i] = p_i(1-q_i) + p_iq_i\mu_i$$
and, because $E[Z_i^2] = \operatorname{Var}(Z_i) + E[Z_i]^2 = \sigma^2 + \mu^2,$
$$E[U_i^2] = \cdots = p_i(1-q_i) + p_iq_i\left(\sigma^2 + \mu^2\right).$$
You can use the equation
$$\operatorname{Var}(W_i) = \operatorname{Var}(aU_i) = a^2 \operatorname{Var}(U_i) = a^2\left(E[U_i^2] - E[U_i]^2\right).$$
At this point it's unclear what you want, because in your code the parameters are themselves random variables.  Are you looking for the unconditional means and variances or, perhaps, for the variances conditional on any specific realization of these parameters, or (based on your bootstrap code) the conditional variance of the mean of the $W_i$?  In any case, you can easily continue this calculation by employing basic properties of expectations and variance.
