So-net無料ブログ作成
検索選択

Xの誤差 補遺 [統計]

生態学会自由集会での発表の補遺。
  1. 回帰分析の場合はXは等間隔であるべきとの指摘を受け、本当のX (x0)をseq(-2, 2, 41)で生成するようにした。
  2. 信頼区間の質問があったので、信頼区間(MCMCではベイズ信用区間)を求めた。


結果

モデル1
plot-21.png

モデル2
plot-22.png

モデル3
plot-23.png

モデル4
plot-24.png

基本的にはx0をrnorm(n, 0, 1)としたときと同様だった。信頼区間(あるいは信用区間)の方は、事前分布を無情報事前分布としたモデル1およびモデル4ではOLSとほぼ同様ということで、久保さんが予想したように結局はOLSと等価なモデルということになっているようだ。

コード。モデル2〜4もこれに準じる。
###
### error in x
###
### x0 <- seq(-2, 2, length = 41)
###


source("~/Documents/RData/JAGS_R.R")
jags <- system("which jags", intern = TRUE)

## model
mod <- "## model 1
var
  N,            # sample size
  x[N],	        # measured x
  y[N],         # measured y
  x0[N],        # real x
  mu.y[N],
  ## parameters
  beta, beta.x,
  tau.x, tau,
  sigma.x, sigma

model {
  for (i in 1:N) {
    y[i] ~ dnorm(mu.y[i], tau);
    mu.y[i] <- beta + beta.x * x0[i];
    x[i] ~ dnorm(x0[i], tau.x);

    ## prior
    x0[i] ~ dnorm(0.0, 1.0E-4);
  }

  ## priors
  beta ~ dnorm(0.0, 1.0E-4);
  beta.x ~ dnorm(0.0, 1.0E-4);
  tau.x ~ dgamma(1.0E-3, 1.0E-3);
  tau ~ dgamma(1.0E-3, 1.0E-3);

  sigma.x <- 1/sqrt(tau.x);
  sigma <- 1/sqrt(tau);
}
"
write(mod, file = model.file)

##
## Monte Carlo
##

set.seed(1234)

test2 <- function(n = 41, beta = 0, beta.x = 0.5,
                  sd.x = 0.1, sd.y = 0.1) {
  x0 <- seq(-2, 2, length = n)
  x <- x0 + rnorm(n, 0, sd.x)
  y <- beta + beta.x * x0 + rnorm(n, 0, sd.y)
  
  init1 <- list(beta = 0, beta.x = 0,
                tau.x = 1, tau.y = 1)
  init2 <- list(beta = -10, beta.x = -10,
                tau.x = 10, tau.y = 10)
  init3 <- list(beta = 10, beta.x = 10,
                tau.x = 0.01, tau.y = 0.01)
  
  post <- jags.run(data = list(N = n, x = x, y = y),
                   inits = list(init1, init2, init3),
                   parameters.to.save = c("beta", "beta.x"),
                   model.file = model.file,
                   data.file = "mcmc1-data.R",
                   init.file = "mcmc1-init.R",
                   cmd.file = "mcmc1.cmd",
                   out = "mcmc1CODA",
                   jags = jags,
                   n.chains = 3,
                   n.iter = 51000,
                   n.burnin = 1000,
                   n.thin = 50)
  post.sum <- summary(post)
  c0 <- lm(y ~ x0)
  c0.est <- c0$coefficients["x0"]
  c0.conf <- confint(c0, "x0")
  c1 <- lm(y ~ x)
  c1.est <- c1$coefficients["x"]
  c1.conf <- confint(c1, "x")
  c2.mean <- post.sum$statistics["beta.x", "Mean"]
  c2.cred <- post.sum$quantiles["beta.x", c(1, 5)]
  c(c0.est, c0.conf, c1.est, c1.conf, c2.mean, c2.cred)
}

mcmc21 <- replicate(2000, test2(n = 41, beta = 0, beta.x = 0.5,
                                sd.x = 0.1, sd.y = 0.1))

save(mcmc21, file = "mcmc21.RData")

nice!(0)  コメント(3)  トラックバック(0) 
共通テーマ:日記・雑感

nice! 0

コメント 3

nakano

hiroki様

拝啓、コメント欄にて失礼いたします。
常々より有益な情報の発信どうもありがとうございます。
2010日本生態学会の自由集会”Xの誤差も統計モデル化”における
”ベイズは~”の投影資料とレジュメに大変興味があるのですが、
もしよろしければ再掲載等、お願いできませんでしょうか? 敬具

nakano
by nakano (2016-01-07 10:51) 

nakano

hiroki様

お手数おかけし申し訳ございません。
どうもありがとうございました。
by nakano (2016-01-10 13:53) 

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

トラックバック 0

この記事のトラックバックURL:
※言及リンクのないトラックバックは受信されません。