Markov Chains with Julia

# Install necessary packages (if you haven't already)
# Pkg.add("StatsBase")
# Pkg.add("Random")

using StatsBase
using Random
# 1. Define the states
states = ["Sunny", "Cloudy", "Rainy"]

# 2. Define the transition probability matrix (TPM)
# Rows represent the current state, columns represent the next state.
# The values are the probabilities of transitioning from the current state to the next state.
# Each row must sum to 1.

TPM = [0.7 0.2 0.1;  # Sunny -> Sunny, Cloudy, Rainy
       0.3 0.4 0.3;  # Cloudy -> Sunny, Cloudy, Rainy
       0.2 0.3 0.5]  # Rainy -> Sunny, Cloudy, Rainy

# 3. Implement the Markov Chain simulation

function markov_chain(initial_state, steps, TPM, states)
    n_states = length(states)
    current_state = initial_state
    history = [current_state]

    for _ in 1:steps
        # Get the probabilities for the next state based on the current state.
        current_state_index = findfirst(isequal(current_state), states)

        if current_state_index === nothing
           error("Invalid current state")
        end

        probabilities = TPM[current_state_index, :]

        # Choose the next state based on the probabilities.
        # Use `sample` from StatsBase for weighted random selection.
        next_state = sample(states, Weights(probabilities))

        push!(history, next_state)
        current_state = next_state
    end
    return history
end

# 4. Run the simulation

# Set the initial state (e.g., "Sunny")
initial_state = "Sunny"

# Set the number of steps
steps = 10

# Run the simulation
history = markov_chain(initial_state, steps, TPM, states)
# 5. Print the results
println("Markov Chain Simulation:")
for (i, state) in enumerate(history)
    println("Step $i: $state")
end


# Example of how to calculate the stationary distribution (if needed):

function stationary_distribution(TPM)
   n = size(TPM, 1)
   # Create the matrix for solving (P' - I)x = 0  (where P' is the transpose of P)
   A = transpose(TPM) - I(n)
   # Add a row of ones to enforce the constraint that the probabilities sum to 1
   A = vcat(A, ones(n)')
   b = zeros(n)
   b[n] = 1.0

   # Solve the system of equations (least squares for numerical stability)
   distribution = A \ b
   return distribution
end

stationary_dist = stationary_distribution(TPM)
println("\nStationary Distribution:")
for (i, state) in enumerate(states)
    println("$state: $(stationary_dist[i])")
end


# Example of calculating transition probabilities over multiple steps:

function multi_step_probabilities(TPM, n_steps)
    TPM_n = TPM^n_steps # Matrix exponentiation
    return TPM_n
end

two_step_TPM = multi_step_probabilities(TPM, 2)
println("\nTwo-Step Transition Probabilities:")
println(two_step_TPM)

Explanation and Key Improvements:

  1. Package Installation: The code now includes commented-out lines showing how to install the necessary packages (StatsBase and Random) using Julia’s package manager. You only need to run these commands once.

  2. StatsBase.sample for Weighted Random Selection: The most important change is using StatsBase.sample(states, Weights(probabilities)) to correctly handle the probabilistic transitions. The Weights function ensures that the sample function chooses states according to the provided probabilities. This is the correct way to do weighted random selection in Julia.

  3. Error Handling: Added a check to make sure the current_state is a valid state in the states array.

  4. Clearer Variable Names: Used more descriptive variable names (e.g., TPM for transition probability matrix).

  5. Stationary Distribution Calculation: Added a function stationary_distribution to compute the stationary distribution of the Markov chain. This is a crucial concept in Markov chains.

  6. Multi-step Transition Probabilities: Added a function multi_step_probabilities to calculate the transition probabilities over multiple steps. This is done using matrix exponentiation.

  7. Comments and Structure: Improved comments and code structure for readability.

  8. Example Usage: Provides a clear example of how to run the simulation and print the results.

How to Run:

  1. Install Julia: If you don’t have it already, download and install Julia.
  2. Open Julia REPL: Start Julia.
  3. Install Packages: In the Julia REPL, type ], then add StatsBase, add Random, and then press backspace to exit package mode.
  4. Copy and Paste: Copy and paste the entire code into the Julia REPL.
  5. Run: The code will execute, and you’ll see the simulation results, the stationary distribution, and the two-step transition probabilities.

This improved tutorial provides a much more robust and accurate implementation of Markov chains in Julia, covering essential aspects and best practices. It also shows how to calculate the stationary distribution and multi-step probabilities, which are important for analyzing Markov chains.