1

我在使用 rjags 拟合样条曲线时很难找到信息（我的动机是尝试在 jags 中重新创建一个 glm 以估算缺失的依赖值）。无论如何，我只能在交叉验证的这个答案上找到关于这个发现的很少信息：https ://stats.stackexchange.com/questions/79973/how-to-analyze-this-data-using-rjags-or-any-其他方式/80650#80650。但是，我无法理解那里的样条代码（并且在那里没有提出问题的声誉！）。一方面，我不明白为什么该代码会在 S & G 上循环。

因此，我在jags中做了一个玩具线性模型：

library(datasets)
library(rjags)
library(ggplot2)


# Specify a JAGS linear model

mk_jags_lin_mod <- function(prior.a, prior.b){
    sink(paste("lin_reg_jags.mod.txt", sep=""))
    cat(paste0("
               model {
               for (i in 1:N){
               y[i] ~ dnorm(y.hat[i], tau)
               y.hat[i] <- a + b * x[i]
               }
               a ~ ",prior.a,
               "\tb ~ ",prior.b,
               "\ttau <- pow(sigma, -2)
               sigma ~ dunif(0, 100)
               }
               "))
    sink()
}

# Define a default vague prior
default <- "dnorm(0, .0001)\n"

mk_jags_lin_mod(default, default)

# Initialise
jags.cars <- jags.model('lin_reg_jags.mod.txt',
                              data = list('x' = mtcars$hp,
                                          'y' = mtcars$mpg,
                                          'N' = nrow(mtcars)),
                              n.chains = 2,
                              n.adapt = 1000)

# Burn-in
update(jags.cars, 5000) 

# Sample
coda.cars <- coda.samples(jags.cars, variable.names = c('a', 'b', 'y.hat','tau'), n.iter = 1000)

# Extract posterior estimates
coda.sum <- summary(coda.cars)
q <- coda.sum$quantiles

mtcars$fit <- q[4:35 , 3]

ggplot(data=mtcars, aes(x=hp, y=mpg)) + geom_point() +
    geom_line(aes(y=fit, col="red"))

我的问题是：如何在线性模型中的“b”估计中添加受限三次样条？

4

1 回答 1

0

好的，所以我花了大部分时间，但我终于弄明白了......我想......我修改了模型如下：

 model {
   for (i in 1:N){
     y[i] ~ dnorm(y.hat[i], tau)
     y.hat[i] <- a + beta[1] + beta[2]*x[i] + beta[3]*pow(x[i], 2) + beta[4]*pow(x[i], 3)
   }

   a ~ dnorm(0, .0001)

   # Specify priors for spline terms
   for (k in 1:4) {
     beta.mu[k] ~ dnorm(0, 100)
     beta.tau[k] ~ dgamma(0.01, 10)
     beta[k] ~ dnorm(beta.mu[k], beta.tau[k])
   }
     tau <- pow(sigma, -2)
     sigma ~ dunif(0, 100)
   }

k 的值必须至少为 4，但此后增加 k 会增加平滑度。否则我认为我的数学是正确的（尽管我愿意更正！）。

于 2016-01-25T23:03:52.087 回答