MarkovChain simulate: Match with python version

src/mc_tools.jl

@@ -295,8 +295,115 @@ function mc_sample_path!(mc::MarkovChain, samples::Array)
  Void
 end
 
+"""
+Simulate time series of state transitions of the Markov chain `mc`.
+
+##### Arguments
+
+- `mc::MarkovChain` : MarkovChain instance.
+- `ts_length::Int` : Length of each simulation.
+- `init::Vector{Int}` : Vector containing initial states.
+- `;num_reps::Union{Int, Void}(nothing)` : Number of repetitions of simulation.
+
+##### Returns
+
+- `X::Array{Int}` : Array containing the sample paths as columns, of shape
+(ts_length, k), where k = length(init) if num_reps=nothing, and k =
+length(init)*num_reps otherwise.
+
+"""
+function simulate(mc::MarkovChain,
+ ts_length::Int,
+ init::Vector{Int};
+ num_reps::Union{Int, Void}=nothing)
+ k = length(init)
+ if num_reps == nothing
+ num_reps = 1
+ end
+ k *= num_reps
+ X = Array(Int, ts_length, k)
+ X[1, :] = repmat(init, num_reps)
+
+ simulate!(mc, X)
+ return X
+end
+
+"""
+Simulate time series of state transitions of the Markov chain `mc`.
+
+##### Arguments
+
+- `mc::MarkovChain` : MarkovChain instance.
+- `ts_length::Int` : Length of each simulation.
+- `init::Int` : Initial state.
+- `;num_reps::Union{Int, Void}(nothing)` : Number of repetitions of simulation.
+
+##### Returns
+
+- `X::Array{Int}` : Array containing the sample path(s) as column(s), of shape
+(ts_length,) if num_reps=nothing; of shape (ts_length, num_reps) otherwise.
+
+"""
+function simulate(mc::MarkovChain,
+ ts_length::Int,
+ init::Int;
+ num_reps::Union{Int, Void}=nothing)
+ init_vec = [init]
+ X = simulate(mc, ts_length, init_vec, num_reps=num_reps)
+ if num_reps == nothing
+ X = vec(X)
+ end
+
+ return X
+end
+
+"""
+Simulate time series of state transitions of the Markov chain `mc`.
+
+##### Arguments
+
+- `mc::MarkovChain` : MarkovChain instance.
+- `ts_length::Int` : Length of each simulation.
+- `;num_reps::Union{Int, Void}(nothing)` : Number of repetitions of simulation.
+
+##### Returns
+
+- `X::Array{Int}` : Array containing the sample path(s) as column(s), of shape
+(ts_length,) if num_reps=nothing; of shape (ts_length, num_reps) otherwise,
+where the initial state is randomly drawn for each simulation.
+
+"""
 function simulate(mc::MarkovChain,
- init::Vector = vcat(1, zeros(n_states(mc)-1)),
- sample_size=1000)
- return mc_sample_path(mc, init, sample_size)
+ ts_length::Int;
+ num_reps::Union{Int, Void}=nothing)
+ if num_reps == nothing
+ init = rand(1:size(mc.p, 1))
+ else
+ init = rand(1:size(mc.p, 1), num_reps)
+ end
+ return simulate(mc, ts_length, init)
+end
+
+"""
+Fill `X` with sample paths of the Markov chain `mc` as columns.
+
+##### Arguments
+
+- `mc::MarkovChain` : MarkovChain instance.
+- `X::Matrix{Int}` : Preallocated matrix of integers to be filled with sample
+paths of the markov chain `mc`. The elements in `X[0, :]` will be used as the
+initial states.
+
+"""
+function simulate!(mc::MarkovChain, X::Matrix{Int})
+ n = size(mc.p, 1)
+ P_dist = [DiscreteRV(vec(mc.p[i, :])) for i in 1:n]
+
+ ts_length, k = size(X)
+
+ for i in 1:k
+ for t in 1:ts_length-1
+ X[t+1, i] = draw(P_dist[X[t]])
+ end
+ end
 end

test/test_mc_tools.jl

@@ -141,6 +141,34 @@ facts("Testing mc_tools.jl") do
  # @fact is_aperiodic(fig) --> true
  end
  end
+
+ context("test simulate shape") do
+ mc = mc3
+ ts_length = 10
+ init = [1, 2]
+ nums_reps = [3, 1]
+
+ @fact size(simulate(mc, ts_length)) --> (ts_length,)
+ @fact size(simulate(mc, ts_length, init)) -->
+ (ts_length, length(init))
+ num_reps = nums_reps[1]
+ @fact size(simulate(mc, ts_length, init, num_reps=num_reps)) -->
+ (ts_length, length(init)*num_reps)
+ for num_reps in nums_reps
+ @fact size(simulate(mc, ts_length, num_reps=num_reps)) -->
+ (ts_length, num_reps)
+ end
+ end
+
+ context("test simulate init array") do
+ mc = mc3
+ ts_length = 10
+ init = [1, 2]
+ num_reps = 3
+
+ X = simulate(mc, ts_length, init, num_reps=num_reps)
+ @fact vec(X[1, :]) --> repmat(init, num_reps)
+ end
 end # facts
 
 end # module