exponentially distributed random number

I want to generate exponentially distributed randm numbers (with mean=lambda) untill the summation of generated numbers is less or equal than T_design=75. thereby, I need to save the random numbers.

[Alex Mcaulley]

See this

Answer 1

[John D'Errico]

Wow. Its been 40 years since I did this sort of stuff. ;-)

I assume you have a stochastic process that you wish to model, as that is what it looks like. So you have some (Poisson) process with exponential inter-arrival time, with rate parameter lambda. You want to generate a set of such arrivals over a span of time T_design.

The simple solution is to just generate those samples, one at a time until you exceed T_design, using a loop.

Be careful here though, because if you use exprnd to generate the samples, exprnd defines the parameter as the mean. Some exponential distribution definitions use the inverse of that. For example, if you read wikipedia, it uses the inverse mean definition of an exponential distribution. (Sigh.)

However, a simple code might look like this:

T_design = 75;
lambda = 3; % I just picked an exponential rate parameter out of a hat
arrivals = [];
while true
  newarrival = exprnd(lambda,1);
  if sum(arrivals) + newarrival < T_design
    arrivals(end+1) = newarrival;
  else
    break
  end
end

So a very simple scheme, that grows a list of arrivals, until you get one that pushes you over the end point.

arrivals
arrivals =
  Columns 1 through 11
       2.6221        9.008      0.30704      0.17039       2.1348       2.1446       3.2566       0.3159       2.9889       6.5892      0.74441
  Columns 12 through 22
       2.8268       4.2603       2.7197       7.0161       6.0755      0.17909      0.13457       1.6591       8.4513       4.3473       3.1225
  Column 23
      0.59099
      
sum(arrivals)
ans =
       71.665

The only flaw to this is it forces you to dynamically grow the array arrivals. That is inefficient. So, perhaps better is to preallocate arrivals to some reasonable size. I expect to see roughly T_design/lambda events over that time span. If I pre-allocate arrivals to be 50% larger than the expected number, I should normally be ok. (We can do better than that estimate, as a function of lambda, of course.) However, this code should be reasonably efficient:

T_design = 75;
lambda = 3; % I just picked an exponential rate parameter out of a hat
arrivals = NaN(1,ceil(1.5*T_design/lambda));
counter = 0;
sumarrivals = 0;
while true
  newarrival = exprnd(lambda,1);
  if sumarrivals + newarrival < T_design
    sumarrivals = sumarrivals + newarrival;
    counter = counter + 1;
    arrivals(counter) = newarrival;
  else
    arrivals = arrivals(1:counter);
    break
  end
end

Can we do better? Of course....

[Mos_bad]

Thank you so much. yeah it is for poisson process. Now I got one more question. I have a while loop whithin a for loop which iterates the lambda values. How I can save the final out put of the nested loops including different generated exponentiall random realization based on different lambda values.

What I am trying to save is "arrivals" for different "m_int" or "lambda_m_int" values.

Here is the code:

M_min=4.5; M_max=8.0; 
m=M_min:0.0001:M_max;
T_design=75; %years
a=4.56; b=1.0; 
alpha=a*log(10);beta=b*log(10);
syms m
F_M=(1-exp(-beta*(m-M_min))) / (1-(exp(-beta*(M_max-M_min)))); % CDF of Mag.
lambda_m=1-F_M; %annual rat
m_int=4.6:0.6:M_max;

for mm=1:length(m_int)
        m_int_=m_int(mm);
        lambda_m_int = double(subs(lambda_m, m,m_int_));
        arrivals = [];
        while true
              newarrival = exprnd(1/lambda_m_int,1);
              if sum(arrivals) + newarrival < T_design
                 arrivals(end+1) = newarrival;
              else
                 break
              end
        end
end

[John D'Errico]

Since the arivals will be a vector of potentially different lengths for every lambda, just save them in a cell array, each set in one cell.