0
$\begingroup$

I have a rv $X$ of the form $$ X=\sum_{i=1}^N Y_i, $$ where $N$ is a discrete rv (often, but not always, Poisson) and $Y_1,\ldots,Y_N$ are continuous random variables, iid and independent from $N$.

I want to compute some quantities with a MC method with a stratified sample. These quantities may depend both on $X$ and on the $X_i$, and I would like to sample $N$ an the $Y_i$ in a way that minimizes the simulation error.

If it is useful to understand the problem, simple examples are: $$ \mathbf E\left[\sum_{i=1}^N \max(0,Y_i-a)\right]\quad\text{or}\quad\mathbf E\left[ \max(0,X-k)\right],$$ more complex ones are of the form $$ \mathbf E\left[\min\left(m,\max\left(0,\sum_{i=1}^N f(Y_1,\ldots,Y_{i-1})\cdot\min(b,\max(0,Y_i-a))\right)-k\right)\right]\quad\text{for some linear $f$.}$$

What I have already tried

  1. Sample $N$ in a stratified way, then sample the appropriate total number of $Y_i$ in a stratified way, shuffle them and attribute them to each simulated scenario. In this case the simulated distribution of $N$ and $Y$ is very close to the theoretical one, but the one for $X$ is still unstable.

  2. Sample $N$ as before, then in separate steps sample, in a stratified way, $Y_1$ for all scenarios where $N\geq1$, then sample $Y_2$ for all scenarios with $N\geq2$, etc. This gives a better approximation for $X$, but a worse one for $Y$.

  3. Sample $N$ as above, Set $\ell$ equal to the max value of $N$, sample $Y_1,\ldots,Y_\ell$ for all scenarios with a latin hypercube sampling method, and then discard the $Y_i$ that I don't need. This method is very slowand similar to 1. in the outcome.

Question

Is there a sampling method that would work for my case?

$\endgroup$
2
$\begingroup$

I'm not really answering your question, but I think you can perform the Latin hypercube sampling quickly, and it does decrease your sampling error for many target metrics. Here is how I would perform the sampling in R:

require(lhs)

set.seed(1980)

N <- 10000
k <- 40 # N and Y_i
X <- randomLHS(N, k)

Y <- X
Y[,1] <- qpois(X[,1], 4)
# check to make sure we don't need a poisson variable bigger than we have
stopifnot(all(Y[,1] < k - 1))

for (i in 2:k)
{
  Y[,i] <- qnorm(X[,i], 2, 1)
}

g <- apply(Y, 1, function(z) 
{
  if (z[1] == 0) return(0)
  return(sum(z[2:(z[1] + 1)]))
})
hist(g, breaks = 25)

# expected mean of g is mean of poisson * mean of iid RV = 4 * 2
mean(g)

################################################################################
# assess sampling error of LHS approach

sims <- 100
res <- numeric(sims)
for (j in 1:sims)
{
  X <- randomLHS(N, k)

  Y <- X
  Y[,1] <- qpois(X[,1], 4)
  # check to make sure we don't need a poisson variable bigger than we have
  stopifnot(all(Y[,1] < k - 1))

  for (i in 2:k)
  {
    Y[,i] <- qnorm(X[,i], 2, 1)
  }
  res[j] <- mean(apply(Y, 1, function(z) 
  {
    if (z[1] == 0) return(0)
    return(sum(z[2:(z[1] + 1)]))
  }))
}
mean(res)
sd(res)

################################################################################
# assess sampling error of SRS approach

res2 <- numeric(sims)
for (j in 1:sims)
{
  X <- rpois(N, 4)
  res2[j] <- mean(sapply(X, function(z) {
    if (z == 0) return(0) else return(sum(rnorm(z, 2, 1)))
  }))
}
mean(res2)
sd(res2)

sd(res2) / sd(res)

```
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.