Poisson distribution simulation  Examples: car accidents in a week, analysis of freeway congestion problem, probability of machine failures in the next week, wait time for a customer at a cache regis

Poisson distribution simulation  Examples: car accidents in a week, analysis of freeway congestion problem, probability of machine failures in the next week, wait time for a customer at a cache register at Trader Joe’s in South Hill, and oodles and noodles more applications.  Assume: The STA bus leaves the Plaza 4 buses an hour on average. But, due to snow, it is totally random that buses arrive at a bus stop. You got up too late and you go to the bus stop as soon as you can, that is, you arrive at the bus stop at a random time. Question: Under this circumstance, what would be your wait time to catch the next bus. A big idea is: bus arrival times 0---b-b---------b---------b---b--b-b-------b---------b--------b---b--599 (in minutes) your arrival time 0------------------------------------y599 your wait time******  generate bunch of bus arrival times and sort small to large simulate your arrival time to the bus stop find your next bus your wait time = your next bus - your arrival time = ****** in the above. Let's formulate STA bus in snow situation. 1.Simulate bus arrival times to the bus stop. 2.Simulate your arrival time to the bus stop. 3.Measure your wait time: next_bus_arrival - your_arrival_time 4.Run the simulation for many times. 5.Find the wait_time distribution.  Now, detail.  Detail 1. Simulate bus arrival times to the bus stop. 1.1We use a random number generator to simulate bus arrival times. ===> the rand( ) function.  The rand( ) function returns a pseudo random number 0 to RAND_MAX (2^31-1 in linux). To generate a random number, rn, bitween 0.0 and 1.0;rn = rand( ) / RAND_MAX.  (by the way, a lot of people do below to create, say, 2 digit random numbers. r_num = rand( ) % 100; since % 100 is 0 to 99. However, this is wrong. The right way of generate 2 digit random number is: divide 0-RAND_MAX in 10 intervals and see where the random number falls. The interval time is, it = RAND_MAX / 100. Then, map it to one of 0 - 99 by the following:  0123.........99 0it2*it3*it99*it to RAND_MAX  If the rand( ) returns a number is between (12*it) and (13*it), the 2 digit random number is 12. )  1.2We will simulate 10 hours (600 minutes). Since four buses are leaving the Plaza every hour on average, forty buses arrives in a ten hour simulation. (If we simulate for only one hour, the distribution of random numbers are nearly uniform.) We will store the arrival time in an array. The unit is measured in minutes:  for(i=0; i< n_buses; i++) r[i] = (float)(599)* ( (float)rand() / (float)RAND_MAX );  That is, 0.0 to 1.0 is mapped to 0.0 to 599 (600 minutes). Thus, r[ ] contains bunch of arrival time measured in minutes. We sort this in small to large order. You use a function to sort:  fsort(&r[0], 600); // or you could do, fsort(r, 600);  Now, r[ ] contains bus-arrival time from the earliest to the latest.  Detail 2. Simulate your arrival time to the bus stop. This is simply;  your_bus_stop_arrival_time = (float)(599)* ( (float)rand() / (float)RAND_MAX );  Detail 3. Measure your wait time Suppose r[ ] contains : 2, 4, 7, 11, 12, 17, 25, …. And, your bus stop arrival time is 8, then your wait time is 11-8 = 3 minutes. You must find your_arrival_time < r[i], then wait_time = r[i] – your_arrival_time.  for( i=0; i<n_buses; i++) if( at <= r[i] ) {wt[j] = (r[i] - your_bus_stop_arrival_time) ; break; }  Detail 4. Run the simulation for many times. What is wt[j] above in Detail 3? We do simulate many times. To do so, we can :  for(j = 0; j< so many times; j++) {Detail 2: your_bus_stop_arrival_time; Detail 3: wt[j] = (r[i] - your_bus_stop_arrival_time) ; }  For each for_loop, we store the wait time in wt[ ].  Detail 5. Find the wait_time distribution. Create a histogram from wt[ ]. That is the distribution of your bus wait time. To do this easier, First sort the wt[ ]. Then, count the number of bus arrival in 5 minute interval. isort(wt, 40); at = 5.0; k=0; for(j=0; j<REPEAT; j++)  { // You arrive at k min at bus stop. // printf("at=%f wt=%f\n", at, wt[j] ); if( (0.0<= wt[j]) && (wt[j] < at) ) hist[0]++; else if( (wt[j] < (2*at)) ) hist[1]++; ……...  else if( (wt[j] < (11*at)) ) hist[10]++; }  OR,  at = 5.0; for(j=0; j<REPEAT; j++) { // You arrive at wt[k] min at bus stop. for(k=0; k<12; k++) { if( (0.0<= wt[j]) && (wt[j] < at) ) { hist[0]++; break; } else if( wt[j] < ((k+1)*at) ) { hist[k]++; break; } } }  hist[k] / total_bus = probability of bus arriving within k*5 min ~ (k+1)*5 min.