1

I am trying to model a system that generates events modeled by a Poisson process.

I am using the following ruby code:

INTERVAL = 0.005
LAMBDA = 167.0
events = Hash.new(0)

def f(x, lambda)
  1 - Math.exp(-lambda * x)
end

random_gen = Random.new
start = Time.now.to_f

while Time.now.to_i - start < 60
  if random_gen.rand < f(INTERVAL, LAMBDA)
    bucket = (Time.now.to_f - start).round.to_i
    events[bucket] += 1
  end

  sleep(INTERVAL)
end

I am trying to generate about 10000 events per minute, or 167 events/second, so I am using $\lambda = 167$.

I am using function f as $1 - e^{(-\lambda x)}$ where $x$ is the interval I am sleeping, I am using this function so that the inter-arrival time of events follows an exponential distribution.

However, I am not getting the expected results, this code generates about 7194 events per minute, with a mean of about 117 events per second. I would expect this code to generate 10000 events per minute with an mean of about 167 events per second.

What am I doing wrong?

Thank you for your help.

UPDATE Fixed typo with sampling time and added random_gen and start definitions

cardinal
  • 24,973
  • 8
  • 94
  • 128
simao
  • 303
  • 3
  • 8
  • I don't understand the point of this implementation. Why are you trying to do this in "real time"? I don't know about the Ruby version of `sleep`, but if it's similar to that in C-like languages, it's both machine dependent and has very low resolution. You can get a "real-time" simulation completely avoiding the sleep call with virtually no change to your code. – cardinal Jul 01 '12 at 23:30
  • I am trying to do this in real time because I need to simulate a distributed system and I need to generate events in real time. How can I get a real time simulation avoiding sleep? – simao Jul 01 '12 at 23:39
  • Just use an extra "counter" variable, say `next_time`, initialized to zero. Then, if `Time.now.to_f - start < next_time`, just continue in the while loop (do nothing). If `Time.now.to_f - start > next_time`, then set `next_time = next_time - log(random_gen.rand)/LAMBDA` and increment your `events` counter. No `sleep` needed! – cardinal Jul 01 '12 at 23:50
  • (Disclaimer: I don't know Ruby, so you may have to adapt the idea to have the right syntax. The `log` in my "code" is the natural logarithm.) – cardinal Jul 01 '12 at 23:52
  • I think you should add your solution as an answer, I think it worked. Still running some tests. – simao Jul 02 '12 at 11:49

3 Answers3

3

How are you getting more than $20$ events per second if you are sleeping for $1/20$ of a second between checks? In case you are actually using intervals of length $0.005$ then you could see up to $200$ events, but you are throwing out each event after the first in each interval. This would lead to an average of $200(1-\exp(-167/200)) = 113.225$ per second, with a binomial distribution instead of a Poisson distribution. You might see slightly fewer if you are actually sampling less often than once every $0.005$ seconds. You should at least test how often your loop is run.

Douglas Zare
  • 10,278
  • 2
  • 38
  • 46
  • 0.05 was a typo yes, I am using 0.005. But why do you say I am throwing out every event after the first? – simao Jul 01 '12 at 22:26
  • I tried running my loop without the if, and indeed about 11000 are generated, meaning that the right condition in the if could fix this – simao Jul 01 '12 at 22:38
  • You are checking whether the count of events in each interval is at least $1$, which occurs with probability $f(,)$. Sometimes a Poisson process will have more than one event on an interval, particularly an interval so long that the average number of events is supposed to be almost $1$, but you are only adding $1$. That's the main reason you have too few. You could reduce this by reducing the length of the interval if you don't want to do anything complicated. A second issue is that it's a bad idea to use sleep. In C# I would make a timer and handle the tick event. – Douglas Zare Jul 01 '12 at 22:58
  • You could have exactly $12000$ checks per minute instead of having about $11000$, sometimes fewer if the system is busy. – Douglas Zare Jul 01 '12 at 22:59
  • Ah I see what you mean. Reducing the interval to 0.001 helps a lot, I got 8105 generated events in one minute, which is closer to what I would expect. Is there any other way to implement this type of simulation? How can I get the number of events that should be generated in the current 0.001 second? Instead of the probability of at least one being generated? – simao Jul 01 '12 at 23:18
  • You can simulate a Poisson distribution by creating a sequence of exponential random variables (the wait times between events), then seeing which partial sum exceeds the length of the interval. If it is the $n$th partial sum the count on that interval is $n-1$. If the length of the interval is small, then on average you won't need to generate many exponential random variables. – Douglas Zare Jul 01 '12 at 23:28
1

You only call f only through f(INTERVAL, LAMBDA) which is a constant equal to 0.9997636. This means that you will go through the if block almost every time.

You sleep a constant time roughly equal to $1/20^{th}$ second so will do a bit less than 20 while cycles per second because of the overhead time of function calls. So you should get around 17-19 events per second.

That answers what is going wrong. Now, to fix it, you need to make f and sleep depend on your random variable. I am not ruby-proficient so this is hard for me to tell, but I suggest something like:

def f()
   # Return an exponential waiting time
   - Math.log(random_gen.rand)
end

and

sleep(f()*0.006)
gui11aume
  • 13,383
  • 2
  • 44
  • 89
  • I am actually using 0.005 for my sleep time The problem with this approach is that sleep doesn't have enough granularity to accept values with that precision, this works for smaller values of $\lambda$, but not for something like $\lambda = 167$ – simao Jul 01 '12 at 22:31
  • 1
    There is no need for the sleep call in the first place and avoiding it will give as near real-time simulation as you're likely to get with Ruby. But, I still don't get why the implementation is done this way. – cardinal Jul 01 '12 at 23:33
0

If you want to generate events from a Poisson process by generating interarrival times you choose X such that $F(X)=U$ where $U$ is uniform on $U[0,1]$ and $F(X)=1-\exp(-\lambda X)$. So you set $U=1-\exp(-\lambda X)$ or $1-U=\exp(-\lambda X)$ or $-\log(1-U)/\lambda =X$. Generate $10000$ $X's$ this way and then count how many $X's$ are less than $60$ seconds. I don't think that is what your code is doing.

Michael R. Chernick
  • 39,640
  • 28
  • 74
  • 143
  • I tried this before, but the problem is that sleep does not have sufficient granularity to support this. – simao Jul 01 '12 at 22:34
  • hm but maybe I can do this every second. Generate 167 Xs and count how many are less than 1 second. Would this work? – simao Jul 01 '12 at 23:25
  • Still doesnt work, I suspect I am doing something wrong when calculating x, it's always much smaller than 1 second, so I always generate 167 events per second. – simao Jul 01 '12 at 23:35