I am trying to generate a large number (approx. 2 millions) of random numbers from an exponential distribution of the form: $$D(t)=\exp(-3000\, t)$$
The sum of the numbers must be 40.
I use the method described in this answer, but I get an underflow error already for 100 numbers. Here is my code:
Select[
Append[#, 40 - Total[#]]& /@
RandomVariate[
ExponentialDistribution[3000],
{1000, 100}
],
Positive[Times @@ #]&
]
Is there a workaround?
@BobHanlon is on the right track with the ErlangDistribution.
First we figure out how to generate a value for $x_1$ given that the sum of the $n$ independent exponential distributed random variables is 40. We have independent random variables $x_1\sim Exponential(\lambda)$ and $\sum_{i=2}^n x_i \sim Erlang(n-1,\lambda)$ and we want to find the conditional distribution of $x_1$ given that $x_1+\sum_{i=2}^n x_i =40$.
The joint distribution is found by
dist = TransformedDistribution[{x1, x1 + x2n}, {x1 \[Distributed] EponentialDistribution[\[Lambda]],
x2n \[Distributed] ErlangDistribution[n - 1, \[Lambda]]}];
The pdf and then cdf of $x_1$ given that the sum is 40 are
pdfx1 = FullSimplify[PDF[dist, {x1, sum}]/Integrate[PDF[dist, {x1, sum}], {x1, 0, sum},
Assumptions -> sum > 0], Assumptions -> n > 1 && x1 > 0 && sum > x1]
(* 40^(1 - n) (-1 + n) (40 - x1)^(-2 + n) *)
cdfx1 = Integrate[pdfx1, {x1, 0, x}, Assumptions -> n > 1 && x > 0 && sum > x]
(* 1 - 40^(1 - n) (40 - x)^(-1 + n) *)
Note that the pdf and cdf do not depend on $\lambda$. We can set the cdf to a Uniform[0,1] random variable and solve for x:
x /. FullSimplify[Solve[cdfx1 == u, x], Assumptions -> 0 <= u <= 1][[1]]
(* 40 - (-40^(-1 + n) (-1 + u))^(1/(-1 + n)) *)
We repeat that process for the remaining values (decrementing the total of 40 as we go).
n = 2000000;
sum = 40;
u = RandomVariate[UniformDistribution[{0, 1}], n - 1];
x = ConstantArray[0, n];
remaining = sum;
Do[x[[i]] = remaining - remaining (1 - u[[i]])^(1/(n - i));
remaining = remaining - x[[i]],
{i, 1, n - 1}];
x[[n]] = sum - Total[x[[1 ;; n - 1]]];
It takes about 7 seconds to generate the 2,000,000 values that sum to 40.
dist = ExponentialDistribution[λ];
Assuming that the random variates are independent, the distribution of the sum of exponential variates is given by the ErlangDistribution
Table[TransformedDistribution @@ {Sum[
x[k], {k, n}], (x[#] \[Distributed] dist) & /@ Range[n]}, {n, 2, 7}]
{ErlangDistribution[2, λ], ErlangDistribution[3, λ],
ErlangDistribution[4, λ], ErlangDistribution[5, λ],
ErlangDistribution[6, λ], ErlangDistribution[7, λ]}
The distribution of the sum of n i.i.d. exponential variates is ErlangDistribution[n, λ]
For n = 2*^6 and λ = 3000
Mean[ErlangDistribution[2*^6, 3000]] // N
(* 666.667 *)
CDF[ErlangDistribution[2*^6, 3000], 40] // N[#, 20] &
(* 4.018024515804428477610^-1627228 )
For the sum to be 40 it is extremely unlikely that the numbers are random
2 million random exponentially distributed numbers with $\lambda=3000$ aren't hard to generate, but their sum is very likely going to be more than 40. So you'll have to settle for fewer numbers - like around 120,000 :
SeedRandom[1];
nums = RandomVariate[ExponentialDistribution[3000], 2*10^6];
pos = First@FirstPosition[Accumulate@nums, x_ /; x > 40]
(* 120026 *)
nums[[1 ;; pos]] // Total
(* 40.0003 *)
Of course, if you sorted them first, you'd get more numbers totaling approximately 40, but then you'd be altering the distribution.
(not a Mathematica answer)
What you are looking for is impossible from a probability point of view. Either you draw a fixed-size i.i.d. sample from some distribution and you have no guarantee on the sum, or you set the size of you sample and its sum but you skew the distribution. Even if you do not set the sample size in advance and generate a sample until the sum has reached a given total, your sample will not correspond to an i.i.d. sample.
Btw, for an exponential distribution of any given parameter, $P(X>40)>0$, therefore even with a sample of size $n=1$ you have to use a truncated distribution.
WorkingPrecision -> 50or similar insideRandomVariate. – MarcoB Oct 09 '21 at 18:22