r

Calculating cumulative precents in R for multiple columns

Suppose I have a discrete time markov chain. I am interested in simulating the Markov Chain over multiple iterations and observing the stationary distribution (i.e. over a long period of time, percent of time spent in all states).

Here is my R code I am currently using.

Here is the Markov Chain:

library(ggplot2)
library(reshape2)  


transition_matrix <- matrix(c(
    0.7, 0.2, 0.1,  # Probabilities of transitioning from A to A, B, C
    0.3, 0.4, 0.3,  # Probabilities of transitioning from B to A, B, C
    0.2, 0.3, 0.5   # Probabilities of transitioning from C to A, B, C
), nrow = 3, byrow = TRUE)


initial_vector <- c(1/3, 1/3, 1/3)

I tried to simulate this as follows:

set.seed(123) 
n_simulations <- 1000
states <- numeric(n_simulations)
current_state <- sample(1:3, 1, prob = initial_vector)

for (i in 1:n_simulations) {
    states[i] <- current_state
    current_state <- sample(1:3, 1, prob = transition_matrix[current_state,])
}

state_names <- c("A", "B", "C")
states_letter <- state_names[states]

df <- data.frame(
    time = 1:n_simulations,
    state = states_letter
)

Finally, I prepare the data for plotting:

cumulative_percentage <- data.frame(
    time = 1:n_simulations,
    A = cumsum(states_letter == "A") / 1:n_simulations * 100,
    B = cumsum(states_letter == "B") / 1:n_simulations * 100,
    C = cumsum(states_letter == "C") / 1:n_simulations * 100
)

cumulative_percentage_melted <- melt(cumulative_percentage, id.vars = "time", 
                                     variable.name = "state", value.name = "percentage")

p2 <- ggplot(cumulative_percentage_melted, aes(x = time, y = percentage, color = state)) +
    geom_line() +
    theme_minimal() +
    labs(title = "Cumulative Percentage of Time Spent in Each State",
         x = "Time Step",
         y = "Cumulative Percentage",
         color = "State") +
    theme(plot.title = element_text(hjust = 0.5)) +
    ylim(0, 100)

p2

state_proportions <- table(states_letter) / n_simulations
print(state_proportions)

enter image description here

Is there a way to change my code so I don't need to manually define A == cumsum, B == cumsum etc etc and have it done for all states?

Solution

  • Your simulation for loop didn't quite make sense to me, so I rolled my own. Probably far from optimal, but the cumulative proportion bit I think is OK, taking advantage of how "vectorised" R is.

    # Transition matrix
states <- LETTERS[1:3]
tran <- matrix(c(
    0.7, 0.2, 0.1,  # Probabilities of transitioning from A to A, B, C
    0.3, 0.4, 0.3,  # Probabilities of transitioning from B to A, B, C
    0.2, 0.3, 0.5   # Probabilities of transitioning from C to A, B, C
), nrow=3, byrow=TRUE, dimnames=list(from=states, to=states))

# Simulation
nsim <- 1e3
sim <- numeric(nsim)
sim[1] <- 1 # Initial state
simseq <- seq_len(nsim)

set.seed(1)
for (i in seq_len(nsim-1)) {
    sim[i+1] <- sample(1:3, 1, prob=tran[sim[i],])
}

# Cumulative proportion
m <- matrix(0, ncol=3, nrow=nsim)
m[cbind(simseq, sim)] <- 1
cumprop <- apply(m, 2, cumsum)/simseq

# Plotting
par(mar=c(3, 3, 1, 1), mgp=c(1.8, 0.6, 0))
matplot(cumprop, type="l", lty=1, xlab="Step n", ylab="Cumulative proportion")
legend("topright", states, lty=1, col=1:3, title="State")

    enter image description here