MC模拟直方图

Question

我正在尝试比较马尔可夫链 (MC) 模拟的直方图和实际数据。我尝试使用下面的代码运行模拟，但我并不完全理解。 R 似乎已经接受了代码，但我不知道如何运行直方图...作为背景，数据是美国经济的扩张和收缩（可在此处找到：http://www.nber.org/cycles.html）。我已经将这两个状态之间的转换矩阵设置为P，其中列总和为 1，并且状态之间的转换计算为“每个状态的转换/月”。我认为n 对应于这里的转换，但我可能是错的......

P <- matrix(c(0.74961, 0.57291, 0.25039, 0.42709),2,2)
P <- t(P)
colSums(P)
n <- 33

MC.sim <- function(n,P) {
    sim<-c()
    m <- ncol(P)    
    sim[1] <- sample(1:m,1) 
    for(i in 2:n){
        newstate <- sample(1:m,1,prob=P[,sim[i-1]])
        sim[i] <- newstate
    }
    sim
}

Answer 1

如 cmets 中所述，您可以使用 hist():

P <- matrix(c(0.74961, 0.57291, 0.25039, 0.42709),2,2)
P <- t(P)
colSums(P)
n <- 33

MC.sim <- function(n,P) {
  sim<-c()
  m <- ncol(P)    
  sim[1] <- sample(1:m,1) 
  for(i in 2:n){
    newstate <- sample(1:m,1,prob=P[,sim[i-1]])
    sim[i] <- newstate
  }
  return(sim)
}

# Save results of simulation
temp <- MC.sim(n,P)

#plot basic histogram
hist(temp)

我没有检查你的模拟。但它只生成 1 或 2 的值。

你也可以使用ggplot():

# Save results of simulation
temp <- MC.sim(n,P)

# Make data frame for ggplot
t <- as.data.frame(temp)
names(t) <- "Sim"

p <- ggplot(t, aes(x = Sim))
p + geom_histogram(binwidth = 0.5)