Tell me more ×
Cross Validated is a question and answer site for statisticians, data analysts, data miners and data visualization experts. It's 100% free, no registration required.

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

share|improve this question
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 1 '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 1 '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 1 '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 1 '12 at 23:52
I think you should add your solution as an answer, I think it worked. Still running some tests. – simao Jul 2 '12 at 11:49

3 Answers

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.

share|improve this answer
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 1 '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 1 '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 1 '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 1 '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 1 '12 at 23:18
show 1 more comment

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.

share|improve this answer
I tried this before, but the problem is that sleep does not have sufficient granularity to support this. – simao Jul 1 '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 1 '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 1 '12 at 23:35

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)
share|improve this answer
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 1 '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 1 '12 at 23:33

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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