2017-10-28 46 views
0

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

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 
} 
+0

您是否在寻找'hist'功能? – akond

+0

该功能对我不起作用,但我可能使用不正确。我只想将仿真绘制为直方图... – Christian

回答

0

正如评论指出,你可以使用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) 
+0

好的,关键的一步似乎是在绘图之前保存仿真。没有意识到我必须这样做。感谢您的帮助。 – Christian

+0

不客气。如果你能将我的答案标记为解决你的问题,那将会很好。谢谢! – FAMG

+0

糟糕!新来Stackoverflow ...答案接受! – Christian