6
$\begingroup$

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?

$\endgroup$
1
  • $\begingroup$ The most immediate source of your problems is a lack of precision to represent the small numbers in your distribution. Add WorkingPrecision -> 50 or similar inside RandomVariate. $\endgroup$
    – MarcoB
    Oct 9, 2021 at 18:22

4 Answers 4

8
$\begingroup$

@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.

$\endgroup$
4
$\begingroup$
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.0180245158044284776*10^-1627228 *)

For the sum to be 40 it is extremely unlikely that the numbers are random

$\endgroup$
2
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

(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.

$\endgroup$
7
  • $\begingroup$ It is not impossible. Conditional distributions are used frequently. $\endgroup$
    – JimB
    Oct 10, 2021 at 3:52
  • $\begingroup$ For example, see stats.stackexchange.com/questions/84098/…. $\endgroup$
    – JimB
    Oct 10, 2021 at 4:26
  • $\begingroup$ But if you constrain the sum then they are not a set of iid random numbers drawn from an exponential distribution. $\endgroup$ Oct 11, 2021 at 5:38
  • $\begingroup$ @AndrewJaffe I interpreted the question as "I am trying to generate a large number (approx. 2 million) of random numbers that sum to 40 and are drawn from the following exponential distribution"; Would this be OK? $\endgroup$ Oct 11, 2021 at 6:48
  • $\begingroup$ @RonaldMonson I think that one way to modify the question in such a way that it does have an answer is suggested by JimB, i.e. how to generate a vector of $n$ numbers from the conditional distribution where the condition is that the sum be 40. Note that the conditional distribution is a distribution of vectors in $\mathcal R^n$ and the vector components are neither iid nor exponentially distributed. $\endgroup$
    – A.G.
    Oct 11, 2021 at 9:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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