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;
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;
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....
Best Answer