0

我使用Stan（特别是 rstan）和贝叶斯单变量线性回归 $y = \beta_0 + \beta_1 x + \varepsilon$。我正在尝试使用 Coda 包来可视化 $\beta$s 的生成轨迹和分布。但是，这会产生错误Error in plot.new() : figure margins too large。 traceplot并且densplot似乎工作正常。问题似乎出在plot.mcmc，它应该产生一个不错的面板输出。您可以在幻灯片“Traceplots and Density Plots”上 看到预期输出的示例。

这是使用数据集的最小（非）工作示例mtcars： 

library(rstan)
library(coda)

stanmodel <- "
data {                     // Data block: Exogenously given information
    int<lower=1> N;        // Sample size
    vector<lower=1>[N] y;  // Response or output. 
                           // [N] means this is a vector of length N
    vector<lower=0, upper=1>[N] x;  // The single regressor; either 0 or 1
}

parameters {               // Parameter block: Unobserved variables to be estimated
    vector[2] beta;        // Regression coefficients
    real<lower=0> sigma;   // Standard deviation of the error term
}

model {                    // Model block: Connects data to parameters
    vector[N] yhat;        // Regression estimate for y
    yhat <- beta[1] + x*beta[2];

    // Priors
    beta ~ normal(0, 10);
    // To plot in R:  plot(function (x) {dnorm(x, 0, 10)}, -30, 30)

    sigma ~ cauchy(0, 5);  // With sigma bounded at 0, this is half-cauchy
                           // http://en.wikipedia.org/wiki/Cauchy_distribution
    // To plot in R:  plot(function (x) {dcauchy(x, 0, 5)}, 0, 10)

    // Likelihood
    y ~ normal(yhat, sigma); // yhat is the estimator, plus the N(0, sigma^2) error
                             // Note that Stan uses standard deviation
}
"

# Designate data
nobs <- nrow(mtcars)
y <- mtcars$mpg
	x <- mtcars$am  # Simple regression version doesn't include constant
data <- list(
    N = nobs,                  # Sample size or number of observations
    y = y,                     # The response or output
    x = x            # The single variable regressor, transmission type
)

# Set a seed for the random number generator
set.seed(123456)

# Run the model
bayes = stan(
    model_code = stanmodel, 
    data = data, # Use the model and data we just defined
    iter = 12000,          # We're going to take 12,000 draws from the posterior, 
    warmup = 2000,         # But throw away the first 2,000
    thin = 10,             # And only keep every tenth draw.
    chains = 3             # But we'll do these 12,000 draws 4 times. 
)

# Use the coda library to visualize parameter trajectories and distributions 
param_samples <- 
    as.data.frame(bayes)[,c('beta[1]', 'beta[2]')]
plot(as.mcmc(param_samples))
4

0 回答 0