2017-09-08 16 views
0

Rの2次元のpfilterのプログラムを作成しました。今度は同じグラフのpfilterとプロットからプロットを描きたいと思います。どうやってするか?同じプロットのデータと一緒にpfilterの結果をプロットする方法

マイコード:

library(magrittr) 
library(pomp) 
library(ggplot2) 
data <- read.csv("C:/Users/admin/Documents/phd/book1.csv", stringsAsFactors = FALSE) 
data <- as.data.frame(data) 

pomp(
    data,times="t", t0=0, 
    rmeasure=Csnippet(" 
        y1 = rnorm(x1+x2,sigma1); 
        y2 = rnorm(x2-x1,sigma2); 
        "), 
    dmeasure=Csnippet(" 
        lik = dnorm(y1,x1+x2,sigma1,1)+dnorm(y2,x2-x1,sigma2,1); 
        lik = (give_log) ? lik : exp(lik); 
        "), 
    rprocess=discrete.time.sim(
    step.fun=Csnippet(" 
         double tx1, tx2; 
         tx1 = rnorm(a11*x1 + a12*x2,nu1); 
         tx2 = rnorm(a21*x1 + a22*x2,nu2); 
         x1 = tx1; x2 = tx2; 
         "), 
    delta.t=1), 
    initializer=Csnippet(" 
         x1 = 0; 
         x2 = 0; 
         "), 
    statenames=c("x1","x2"), 
    paramnames=c("a11","a12","a21","a22","sigma1","sigma2","nu1","nu2"), 
    params=c(a11=0.5,a12=-0.1,a21=0.2,a22=-1,nu1=0.3,nu2=0.1,sigma1=0.1,sigma2=0.3) 
) -> parus 

pf <- pfilter(parus, Np=1000, 
       params=c(a11=0.5, a12=-0.1, a21=0.2, a22=-1, nu1=0.3, nu2=0.1, 
         sigma1=0.1, sigma2=0.3)) 
logLik(pf) 

ggplot(data, aes(x= t, y= y1)) + 
    geom_point() + 
    geom_smooth(method="pfilter", se=FALSE, colour = "red") 

マイデータサンプル:

data <- 
structure(list(t = 1:10, y1 = c(-0.058886245, -0.045854389, 0, 
-0.06009386, 0.053712929, 0, 0.047628037, 0.061425641, 0, 0), 
    y2 = c(0.018451378, -0.012812011, -0.025850852, 0.021278254, 
    -0.008261274, 0.025270473, -0.008357995, -0.0146046, -0.03778035, 
    -0.004632415)), .Names = c("t", "y1", "y2"), class = "data.frame", row.names = c(NA, 
-10L)) 
+0

あなたの例が使えるようにいくつかのサンプルデータを 'dput'してください。 –

答えて

0

あなたはPFからデータを抽出し、ggplotにそれらを追加構築する必要があります。例えば :あなたが欲しいことがあり、残りの行を追加することができます。同様に

pfdata<- pf$data %>% as.data.table(keep.rownames = T) %>% melt(idvars = rn) 
pfdata[rn == "y1",t:=1:10 ] 
pfdata[rn == "y2",t:=1:10 ] 

pfdata <- pfdata %>% dcast.data.table(variable + t ~rn) 


data <- merge(data, pfdata, by = "t") 

ggplot(data) + 
    geom_point(aes(x= t, y= y1.x)) + 
    geom_line(aes(x= t, y= y1.y))+ 
    geom_line(aes(x= t, y= y2.y)) 

関連する問題