We have a few parameters to work with:
delay)shift)shape and rate)We can use the rgamma function in R to generate a random sample from the Gamma distribution with the shape and rate parameters.
n = 10000
shape = 1
rate = 2
shift = 0.5
interboxinterval <- rgamma(n,shape,rate)+shift
hist(interboxinterval,breaks=seq(0,7,by=0.5),main="Time between boxes",xlab="Inter-box Interval")
We’ll define a function that simulates the number of boxes we will see in myTime seconds. Because we don’t know about the gamma function, we’re not able to simulate from within R.
Frustratingly, we are going to resample from the sample given below, which was generated by Javascript.
myrand =
c(0.542602579928852,
1.0538298255630176,
0.6363779415771225,
1.688190677252905,
1.6783066727594513,
2.6893253731961573,
0.7888557981435277,
1.6104436295012419,
0.17256034313243826,
0.4558008915297854,
0.2032839729564929,
0.639712686741254,
0.7036194213288965,
0.2936640641229919,
1.4907274543772233,
0.05063973852869065,
1.5888769898984827,
0.4858788492377629,
0.9524655134110718,
0.8820206766687115,
0.8240021317282659,
1.1674574327966813,
0.7228057828145512,
0.7799503979300636,
1.4628091543927393,
0.40961087708610383,
1.4353589514226104,
0.5987368944934536,
2.211404379334857,
1.8898642568932562,
0.2708042525020506,
0.09770222195104074,
0.0739909105792561,
0.8295848027015739,
0.35310268174701703,
0.5089461709410353,
0.020796093794425676,
0.6125450879807594,
0.5336101931513086,
0.49992206015359697,
1.7503044333244797,
0.12434279640253881,
0.6033522259128451,
0.589929359002586,
0.274633060890688,
2.1195930113469825,
1.3102330192654474,
0.598367753165459,
4.435521345898014,
0.636901699176378,
0.954893143239095,
1.3646683313090684,
0.1476180745106118,
0.015605854808236225,
0.4390463315816505,
0.9164478992769152,
0.37258603493420844,
0.12158819231752493,
0.7082855189364727,
1.8026465026321064,
1.5233530576295862,
0.45747216833189636,
0.02018198774188074,
2.521523526138969,
1.807386965807478,
0.08609240681746673,
0.9198013694916899,
1.1467167156344833,
0.5464402105184519,
1.7110069060587652,
1.4664125073773762,
2.198088370936873,
0.5923942808630472,
2.059607535252995,
0.7233365195212779,
0.12617175225144944,
0.3445536395192657,
4.818531818938236,
0.8666071510613672,
0.27335991797896686,
3.5842389576749563,
0.41124728352951145,
1.1163628215575871,
0.5753220895203456,
0.34142670050451823,
1.407568642625773,
1.1977206574333281,
0.8373914521970929,
1.2285192161481102,
1.6794689810642303,
0.6731812192833158,
0.9452466124147185,
4.960916431538967,
0.19005538072681555,
0.7418607258758575,
0.617034041137649,
0.9478174024574951,
0.7172509196961968,
1.745623817299884,
0.014523426925867347,
0.02302214598829502,
0.7654743155712774,
2.259452450936171,
0.0804732669560493,
1.3734533597624494,
0.5366073406254775,
1.5465582463943048,
0.9886044675657534,
2.252959211367411,
3.36076569026019,
0.18522809442875052,
0.47724764294294314,
0.7678175057063028,
0.2712040775196358,
2.1236135548067936,
0.0641498206202489,
1.4445815906374595,
0.022238692055816835,
5.357634528279393,
0.6898148458375476,
2.0307374140214622,
0.017881557253850624,
0.4587582566864349,
1.515793365524987,
0.04920461300887018,
1.8350012762729442,
0.7980065511747395,
3.5510744762869852,
2.9790493493937977,
2.188295788641627,
3.4162369666283685,
4.833698731169567,
0.4264713591311623,
2.80629028466563,
1.6979710440008953,
0.41703982896883934,
0.21108059750190572,
6.092520418878765,
0.43458952913661575,
2.6675727146310337,
0.12665407493815486,
1.0184499847720407,
0.4298657553164043,
0.19417528364802533,
0.46700876329669333,
1.2801835865014173,
2.499592238733759,
0.35978843062672305,
0.4259265655028755,
0.7694869024795504,
0.24444578754464133,
0.30113177608272196,
0.11719685901559274,
0.3076390828001878,
1.9329476853531324,
0.17067635051938654,
1.919761452293006,
1.5475999794013575,
2.075742181100213,
0.3850213240976466,
2.63023023555751,
4.188541757132042,
2.6656058682837096,
0.4746563755162958,
0.4687500646742593,
1.9774854028805922,
1.0825268724623274,
1.519965491414461,
3.086489606573306,
0.17785071292083474,
1.519160756811418,
0.41986850558618205,
0.6981456112796207,
1.8674641889650925,
0.5272490637655568,
1.3449637568087556,
1.4280299446322675,
0.2733123663339775,
0.280804897164523,
0.7737442376302235,
1.194728359745683,
0.30179113682161035,
2.856309317456353,
1.1943961038135302,
0.6498420378048678,
0.2648314824201625,
0.1993727165082136,
0.3557854211851964,
0.39463998319865173,
4.487181331602883,
0.5433629848823083,
1.642459075887904,
0.45591705176432507,
2.4446165064128436,
0.8159744208825696,
0.7119394647123211,
2.4590540079488554,
0.284593930926141,
0.047239281915853255)
boxesUntilT <- function (myTime,shape,rate,shift,delay) {
n = 500 # the number of draws from the distribution
#nsecs <- rgamma(n,shape,rate)+shift
nsecs <- sample(myrand,n,replace=T)+shift
boxsecs <- cumsum(nsecs+delay)
w = which(boxsecs<myTime)
w[length(w)]
}
Then, we can run this experiment 1000 times to get an estimate of how many boxes we will see on average in, say, 120 seconds.
shape = 1
rate = 2
shift = 0.5
delay = 1
myTime = 120
results <- replicate(1000,boxesUntilT(myTime,shape,rate,shift,delay))
mean(results)
## [1] 44.82
hist(results,main="Histogram of boxes seen at 120 Seconds",xlab="Number of Boxes")
We can test a range of parameters using a for loop. For example, here we try values of shift of 0.1, 0.15, 0.2, 0.25,…0.5.
df <- data.frame(time = numeric(0),boxes=numeric(0))
shape = 1
rate = 2
shift = 0.5
delay = 1
myTime = 30
for (shift in seq(0.1,0.5,0.05)) {
df <- rbind(df,data.frame(shift=shift, boxes = mean(replicate(100,boxesUntilT(myTime,shape,rate,shift,delay)))))
}
print(df)
## shift boxes
## 1 0.10 12.89
## 2 0.15 12.35
## 3 0.20 12.49
## 4 0.25 11.95
## 5 0.30 11.75
## 6 0.35 11.63
## 7 0.40 11.36
## 8 0.45 11.00
## 9 0.50 10.75
We can even try sweeping across multiple parameters at once.
df <- data.frame(time = numeric(0), rate=numeric(0),boxes=numeric(0))
shape = 1
rate = 2
shift = 0.5
delay = 1
myTime = 30
for (shift in seq(0.1,0.5,0.05)) {
for (rate in seq(1,3,0.5)) {
df <- rbind(df,data.frame(shift=shift, rate=rate,boxes = mean(replicate(100,boxesUntilT(myTime,shape,rate,shift,delay)))))
}
}
print(df)
## shift rate boxes
## 1 0.10 1.0 12.78
## 2 0.10 1.5 13.07
## 3 0.10 2.0 12.78
## 4 0.10 2.5 12.65
## 5 0.10 3.0 12.81
## 6 0.15 1.0 12.61
## 7 0.15 1.5 12.91
## 8 0.15 2.0 12.64
## 9 0.15 2.5 12.76
## 10 0.15 3.0 12.41
## 11 0.20 1.0 12.51
## 12 0.20 1.5 12.44
## 13 0.20 2.0 12.33
## 14 0.20 2.5 12.22
## 15 0.20 3.0 12.37
## 16 0.25 1.0 11.98
## 17 0.25 1.5 11.98
## 18 0.25 2.0 12.21
## 19 0.25 2.5 11.89
## 20 0.25 3.0 12.57
## 21 0.30 1.0 11.66
## 22 0.30 1.5 12.00
## 23 0.30 2.0 11.35
## 24 0.30 2.5 12.01
## 25 0.30 3.0 11.91
## 26 0.35 1.0 11.52
## 27 0.35 1.5 11.10
## 28 0.35 2.0 11.40
## 29 0.35 2.5 11.31
## 30 0.35 3.0 11.88
## 31 0.40 1.0 11.00
## 32 0.40 1.5 11.03
## 33 0.40 2.0 11.18
## 34 0.40 2.5 11.33
## 35 0.40 3.0 11.14
## 36 0.45 1.0 11.23
## 37 0.45 1.5 11.18
## 38 0.45 2.0 11.00
## 39 0.45 2.5 11.27
## 40 0.45 3.0 11.06
## 41 0.50 1.0 10.95
## 42 0.50 1.5 10.91
## 43 0.50 2.0 10.90
## 44 0.50 2.5 10.98
## 45 0.50 3.0 10.85