主页 / coding / 问题 / 79075411
Accepted
farrow90
Asked: 2024-10-11 00:45:09 +0800 CST

计算 R 中多列的累计百分比

  • 772

假设我有一个离散时间马尔可夫链。我感兴趣的是模拟马尔可夫链的多次迭代并观察其平稳分布（即在一段较长的时间内，处于所有状态的时间百分比）。

这是我目前正在使用的 R 代码。

这是马尔可夫链：

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)

我尝试按如下方式模拟此过程：

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
)

最后，我准备绘图的数据：

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)

在此处输入图片描述

有没有办法改变我的代码，这样我就不需要手动定义 A == cumsum、B == cumsum 等等，并对所有状态完成它？

1 个回答

  1. Best Answer

    您的模拟for循环对我来说不太合理，所以我自己做了一个。可能远非最优，但我认为累积比例位还可以，充分利用了 R 的“矢量化”方式。

    # 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")

    在此处输入图片描述

    • 2

相关问题