Re: Discrepancy Between 2 Different Ways Of Calculating - What's Wrong, Please?
- From: "Daniel Lichtblau" <danl@xxxxxxxxxxx>
- Date: 28 Apr 2006 12:33:54 -0700
Mark Lawton wrote:
I am trying to answer the following question: if an event is expected
to happen n times (e.g. twice) in a 3 year period, what is the
probability of it happening 4 or more times?
Here's a Poisson distribution function I wrote:
(%i46) poisso(x,lam) := (i: x, sum((%e^(lam * -1) * lam^i)/i!, i,0,i));
lam (- 1) i
%e lam
(%o46) poisso(x, lam) := (i : x, sum(----------------, i, 0, i))
i!
(%i40) poisso(4,2),numer;
(%o40) 0.94734698265629
So, there's about a 5% chance of it happening 4 times. This can be
checked with the Poisson Distribution calculator at
http://www.ciphersbyritter.com/JAVASCRP/BINOMPOI.HTM .
I also wanted to be able to calculate it day by day. Assuming that 3
years = 1095 days, if the event is expected to happen twice, then its
daily probability will be 2/1095, or 0.001826. Using the binomial
calculator at http://www.ciphersbyritter.com/JAVASCRP/BINOMPOI.HTM ,
with p=0.001826, n=1095 and k=4, we again get an answer of about 5%.
I then wrote a monte-carlo simulator (below) - and it consistently
gives answer of around 14-15%. I cannot see where the discrepancy is
arising from - can anyone spot it for me, please?
In this Maxima program, 1000 simulations of 3 year periods (1095 days)
are done. A day is "counted" if a random integer in the range 0..1095
equals 0 or 1 (which leads to an expected value of 2 counts in the 1095
days), and an "event" is counted if 4 or more counts occur in the 1095
days.
(events: 0,
for outer: 1 thru 1000 do
(count: 0,
for inner: 1 thru 1095 do
if random(1095) <= 1 then
count: count + 1,
if count >= 4 then
events: events + 1),
print("Event happened 4 or more times in", float(events/10), "% of the
tests"))$
(I also wrote it in visual basic, and got the same answer):
Dim count, events As Integer
Dim level As Double
level = 2/1095
events = 0
For outer = 1 To 1000
count = 0
For inner = 1 To 1095
If Rnd() < level Then
count = count + 1
End If
Next
If count >= 4 Then
events = events + 1
End If
Next
Messagebox Cstr(events/10) & "%, (" & Cstr(events) & ")"
You are working with two different distributions. In your simulation
you assume, in effect, that the events are independent and (discrete)
uniformly distributed. With this distribution the probability of four
or more occurrences could be found as below (using Mathematica).
n = 1095;
p = 2/n;
ss = Sum[Binomial[n,j]*p^j*(1-p)^(n-j), {j,4,n}];
In[8]:= N[ss]
Out[8]= 0.142711
This is in accord with your simulation. And mine (below).
In[4]:= runs1k = Table[Apply[Plus,UnitStep[
-Table[Random[Integer,{0,1095}],{1095}]+3/2]], {1000}];
In[5]:= Length[Select[runs1k, #>=4&]]
Out[5]= 135
In[6]:= runs1k = Table[Apply[Plus,UnitStep[
-Table[Random[Integer,{0,1095}],{1095}]+3/2]], {1000}];
In[7]:= Length[Select[runs1k, #>=4&]]
Out[7]= 156
In[8]:= runs1k = Table[Apply[Plus,UnitStep[
-Table[Random[Integer,{0,1095}],{1095}]+3/2]], {1000}];
In[9]:= Length[Select[runs1k, #>=4&]]
Out[9]= 143
Daniel Lichtblau
Wolfram Research
.
- References:
- Prev by Date: Discrepancy Between 2 Different Ways Of Calculating - What's Wrong, Please?
- Next by Date: Dream ƒrequency Mathematics Applications
- Previous by thread: Discrepancy Between 2 Different Ways Of Calculating - What's Wrong, Please?
- Next by thread: Re: Discrepancy Between 2 Different Ways Of Calculating - What's Wrong, Please?
- Index(es):
Relevant Pages
|
|
