誤差のある説明変数

Kubolog 2009-10-12にあった話。自分でもやってみる。

set.seed(1234)

test1 <- function(beta = 0, beta.x = 0.5,
			n = 20,
			sd.x = 0.1, sd.y = 0.1) {
  x0 <- rnorm(n, 0, 1)
  x <- x0 + rnorm(n, 0, sd.x)
  y <- beta + beta.x * x0 + rnorm(n, 0, sd.y)
  c0 <- lm(y ~ x0)$coefficients[2]
  c1 <- lm(y ~ x)$coefficients[2]
  c(c0, c1)
}

t1 <- replicate(5000, test1(n = 100))

pdf("Rplot.pdf", width=512/72, height=512/72)
plot(density(t1[1,]), col = "black", las = 1, main = "density")
lines(density(t1[2,]), col = "red")
dev.off()

結果

と、傾きが小さい方に偏る。

では、ベイズ推定ならうまくいくかというと

library(rjags)

set.seed(1234)
n <- 100	# sample size
sd.x <- 0.1
sd.y <- 0.1
b.x <- 0.5
x0 <- rnorm(n, 0, 1)
x <- x0 + rnorm(n, 0, sd.x)
y <- b.x * x0 + rnorm(n, 0, sd.y)

model <- jags.model("x_error.bug",
                    data = list(N = n,
                                x = x, y = y),
                    n.chains = 3, n.adapt = 10000)
post <- coda.samples(model,
                     variable = c("beta", "beta.x",
                                  "sigma.x", "sigma.y"),
                     n.iter = 20000, thin = 20)

pdf("Rplot2.pdf", width=512/72, height=768/72)
plot(post)
dev.off()

モデル

var
  N,  # sample size
  beta, beta.x,
  tau.x, tau.y,
  sigma.x, sigma.y, 
  x[N], y[N], x0[N], y0[N]

model {
  for (i in 1:N) {
    y[i] ~ dnorm(y0[i], tau.y);
    y0[i] <- beta + beta.x * x0[i];
    x[i] ~ dnorm(x0[i], tau.x);
    ## prior
    x0[i] ~ dnorm(0.0, 1.0E-6);
  }

  beta ~ dnorm(0.0, 1.0E-6);
  beta.x ~ dnorm(0.0, 1.0E-6);
  tau.x ~ dgamma(1.0E-4, 1.0E-4);
  tau.y ~ dgamma(1.0E-4, 1.0E-4);
  sigma.x <- 1/sqrt(tau.x);
  sigma.y <- 1/sqrt(tau.y);
}

結果

sigma.xとsigma.yが不安定。X方向で調整するか、Y方向で調整するかで局所的に安定になるのであろう。ということで、このモデルではあまりうまくいかない。

タグ：MCMC