Combing $n$ 2-point distributions I have $X_1,X_2,\dots,X_{40}$ independent distributions of the following form:
$f_X(x) = \begin{cases}
0 &\text{with probability }p=0.4\\
10 &\text{with probability }q=0.6
\end{cases}$
I also have $Y_1,Y_2,\dots,Y_{10}$ independent distributions of the following form:
$f_Y(y) = \begin{cases}
0 &\text{with probability }p=0.94\\
100 &\text{with probability }p=0.06
\end{cases}$
Essentially, I'd like to find a distribution for $X_1 + \dots + X_{40} + Y_1 + \dots + Y_{10}$. My approach right now is to scale each by $1/10$ and $1/100$ respectively to turn this into a series of independent Bernoulli distributions, so that in the end I'm summing two binomial distributions, and can use convolution.
However, I've realized that the problem then becomes not knowing how to "rescale" as it were. I'm wondering if anybody can think of a good way to go about doing this. Ultimately, I'd like to use these distributions in a Monte Carlo simulation.
 A: Let's just look at the $X$s for a start:
As I understand your question, $X_i \sim 10\mathsf{Bern}(.6).$ So
$S = \sum_{i=1}^{40} X_i = 10\mathsf{Binom}(40, .6).$
Then $E(S) = 10np = 10[40(.6)] = 240.$ Also,
$Var(S) = 10^2np(1-p) = 100(9.6) = 960$ and $SD(S) \approx 30.98.$
You mention simulation. In R statistical software one method of
simulating $S$ is shown below. A sample of million realizations of $S$ should
have sample mean and SD approximating the numbers above to roughly three
digit accuracy. (Of course every run of the simulation will yield slightly
different results.)
s = 10*rbinom(10^6, 40, .6)
mean(s);  sd(s)
[1] 239.9632  # aprx E(S) = 240
[1] 30.97314  # aprx SD(S) = 30.98

Of course, $S$ is discrete, taking values that are multiples of $10,$ but it should be approximately normal.
Another way to simulate a million realizations of $S$ is to use the
sample function in R to sample directly from integers $0$ and $10$ with the specified probabilities.
Results are about the same:
s = replicate(10^6, sum(sample(c(0,10), 40, repl=T, prob=c(.4,.6))))
mean(s);  sd(s)
[1] 240.0496
[1] 30.99703

In the histogram below, each achieved value of $s$ is located at the
center of a histogram bar; the density curve is for $\mathsf{Norm}(\mu-240, \sigma = 30.98).$

If you do the theory for the $Y$s, you can get the mean and standard
deviation of $T = \sum_{i=1}^{10} Y_i.$ Analogous simulations should work OK.
Because this distribution sums fewer Bernoulli random variables and has
wider gaps between values, you will have to check whether a normal
approximation seems serviceable for your purposes. If so, then $S + T$ (a sum of multiples of binomials) should also
be "approximately" normal.
