using SimJulia
using Distributions
using RandomStreams
using Distributions
const SEED = 12345
rand_dist(rng::MRG32k3a, Dist::Distribution) = quantile(Dist, rand(rng))
seeds = [SEED, SEED, SEED, SEED, SEED, SEED]
gen = MRG32k3aGen(seeds)
include("tally.jl")
include("tallystore.jl")
println("M/M/1 with processes")
#rates
λ = 1.9
μ = 2.0
ρ = λ/μ
type System
W
arrival
service
counter
arrgen
servgen
function System(text::String)
s = new()
s.W = Tally(text)
s.arrival = Exponential(1.0/λ)
s.service = Exponential(1.0/μ)
s.arrgen = next_stream(gen)
s.servgen = next_stream(gen)
return s;
end
end
function restart(s::System)
init(s.W)
next_substream!(s.arrgen)
next_substream!(s.servgen)
end
s = System("Waiting Times")
# Allow a simulation with a fixed number of clients
function source(env::Environment, s::System, counter::Resource, limit::Float64)
i = 0
while (true)
yield(Timeout(env, rand_dist(s.arrgen, s.arrival)))
if (now(env) > limit) break end
i += 1
Process(env, customer, s, i, counter)
end
end
function source_fixed(env::Environment, s::System, counter::Resource, nCusts::Int64)
for i = 1:nCusts
yield(Timeout(env, rand_dist(s.arrgen, s.arrival)))
Process(env, customer, s, i, counter)
end
end
function customer(env::Environment, s::System, idx::Int, counter::Resource)
# Record the arrival time in the system
arrive = now(env)
yield(Request(counter))
# The simulation clock now contains the time when the client goes to the server.
wait = now(env) - arrive
# Record the waiting time
add(s.W, wait)
yield(Timeout(env, rand_dist(s.servgen, s.service)))
yield(Release(counter))
end
function onesim(s::System)
env = Environment()
s.counter = Resource(env, 1)
Process(env, source_fixed, s, s.counter, 10000)
run(env)
end
onesim(s)
average(s.W)
variance(s.W)
meanWaits = TallyStore("Temps d'attente moyen")
Nobs = 100
for n=1:Nobs
restart(s)
onesim(s)
w = average(s.W)
N = nobs(s.W)
add(meanWaits, w)
println("Sim $n. Average waiting time: $w. Number of observations: $N")
end
average(meanWaits.t)
variance(meanWaits.t)
stdev(meanWaits.t)
var(meanWaits.obs)
function ci_normal(n::Int64, mean::Float64, stdev::Float64, α::Float64)
z = quantile(Normal(), 1-α/2)
w = z*stdev/sqrt(n)
# Lower bound
l = mean - w
# Upper bound
u = mean + w
return l, u
end
ci_normal(nobs(meanWaits.t), average(meanWaits.t), stdev(meanWaits.t), 0.05)
function ci_tdist(n::Int64, mean::Float64, stdev::Float64, α::Float64)
z = quantile(TDist(n-1), 1-α/2)
w = z*stdev/sqrt(n)
# Lower bound
l = mean - w
# Upper bound
u = mean + w
return l, u
end
ci_tdist(nobs(meanWaits.t)-1, average(meanWaits.t), stdev(meanWaits.t), 0.05)
Average waiting time according to M/M/1 formulas. In stationnary regime, the mean waiting time is $$ \overline{W} = \frac{\frac{\rho}{\mu}}{1-\rho} = \frac{\rho}{\mu-\lambda}. $$
w = (ρ/μ)/(1-ρ)
sortedObs = sort(meanWaits.obs)
Bootstrap?
# Construct a bootstrap sample
unif=next_stream(gen)
m = 200
wboostrap = TallyStore("Boostrap estimator")
for j = 1:m
xb = Array(Float64,Nobs)
for i=1:Nobs
k = Int64(ceil(rand(unif)*Nobs))
xb[i] = meanWaits.obs[k]
end
add(wboostrap, mean(xb))
end
average(wboostrap.t)
sortwb = sort(wboostrap.obs)
stdev(meanWaits.t)
α = 0.05
ceil(m*(1-α/2))
ceil(m*α/2)
y = average(meanWaits.t)
2*y-sortwb[195]
2*y-sortwb[5]