Lindley recurrence

Add the unofficial RandomStreams package if needed

In [1]:
Pkg.clone("https://github.com/prsteele/RandomStreams.jl")
Pkg.update()

INFO: Cloning RandomStreams from https://github.com/prsteele/RandomStreams.jl
INFO: Computing changes...
INFO: Updating METADATA...
INFO: Updating RandomStreams master...
INFO: Computing changes...
INFO: No packages to install, update or remove
In [1]:
using Distributions
using RandomStreams  # provides the MRG32K3a generator along with random streams
In [2]:
Pkg.update()

INFO: Updating METADATA...
INFO: Updating cache of DataFrames...
INFO: Updating cache of Rmath...
INFO: Updating cache of DataArrays...
INFO: Updating cache of Conda...
INFO: Updating cache of Compat...
INFO: Updating cache of PDMats...
INFO: Updating cache of DataFrames...
WARNING: Package RandomStreams: skipping update (dirty)...
INFO: Computing changes...
INFO: Upgrading Compat: v0.9.2 => v0.9.3
INFO: Upgrading Conda: v0.3.2 => v0.4.0
INFO: Upgrading DataArrays: v0.3.9 => v0.3.10
INFO: Upgrading PDMats: v0.5.0 => v0.5.1
INFO: Upgrading Rmath: v0.1.3 => v0.1.4
INFO: Building Conda
INFO: Building Rmath
In [3]:
const SEED = 12345

seeds = Array[[SEED, SEED, SEED, SEED, SEED, SEED]]
gen = MRG32k3aGen(seeds[1])
arrival_gen = next_stream(gen)
Out[3]:
Full state of MRG32k3a generator:
Cg = [12345,12345,12345,12345,12345,12345]
Bg = [12345,12345,12345,12345,12345,12345]
Ig = [12345,12345,12345,12345,12345,12345]
  • Cg is the current state of the generator
  • Bg is the first state of the substream
  • Ig is the first state of the stream
In [4]:
service_gen = next_stream(gen)
Out[4]:
Full state of MRG32k3a generator:
Cg = [3692455944,1366884236,2968912127,335948734,4161675175,475798818]
Bg = [3692455944,1366884236,2968912127,335948734,4161675175,475798818]
Ig = [3692455944,1366884236,2968912127,335948734,4161675175,475798818]
In [5]:
rand_dist(rng::MRG32k3a, Dist::Distribution) = quantile(Dist, rand(rng))
Out[5]:
rand_dist (generic function with 1 method)

Lindley recurrence: $$ W_1 = 0,\quad W_{i+1} = \max(0,\; W_i + S_i - A_{i+1}). $$

In [7]:
function mean_wait()
    N = 1000000 # Number of clients
    Waits = Array(Float64, N)
    Waits[1] = 0.0 # The first client does not wait.
    lambda = 1.9   # arrival rate
    mu = 2.0       # service rate
    arrival = Exponential(1.0/lambda)
    service = Exponential(1.0/mu)

    for i in 2:N
        Waits[i] = max(0, Waits[i-1]-rand_dist(arrival_gen, arrival)+rand_dist(service_gen, service))
    end
    
    return Waits
end

WARNING: Method definition mean_wait() in module Main at In[6]:2 overwritten at In[7]:2.
Out[7]:
mean_wait (generic function with 1 method)
In [8]:
Waits = mean_wait();
In [9]:
mean(Waits)
Out[9]:
9.957419540675662
In [10]:
next_substream!(arrival_gen)
In [11]:
arrival_gen
Out[11]:
Full state of MRG32k3a generator:
Cg = [870504860,2641697727,884013853,339352413,2374306706,3651603887]
Bg = [870504860,2641697727,884013853,339352413,2374306706,3651603887]
Ig = [12345,12345,12345,12345,12345,12345]
In [12]:
next_substream!(service_gen)
In [13]:
service_gen
Out[13]:
Full state of MRG32k3a generator:
Cg = [3119395571,2178405402,1065030501,3980307777,2117495919,1836828492]
Bg = [3119395571,2178405402,1065030501,3980307777,2117495919,1836828492]
Ig = [3692455944,1366884236,2968912127,335948734,4161675175,475798818]
In [15]:
Waits = mean_wait()
mean(Waits)
Out[15]:
0.4487765693750297
In [16]:
arrival_gen
Out[16]:
Full state of MRG32k3a generator:
Cg = [2826986986,2700790508,2240279461,1099791616,3966731335,486127084]
Bg = [870504860,2641697727,884013853,339352413,2374306706,3651603887]
Ig = [12345,12345,12345,12345,12345,12345]
In [14]:
Waits
Out[14]:
100-element Array{Float64,1}:
 0.0      
 1.17113  
 0.828511 
 0.0      
 0.694298 
 0.653771 
 0.150675 
 0.0      
 0.0      
 0.0      
 0.533166 
 0.0      
 0.265602 
 ⋮        
 0.0      
 0.0      
 0.238497 
 0.0530849
 0.0      
 0.0      
 0.556536 
 0.0      
 0.160575 
 0.0      
 0.0      
 0.0
In [ ]: