r - 使用 Quadprog 对系数范围/限制的线性模型的约束

Question

我正在做一些实验，我当然知道为什么很少需要约束系数，但是这里有。

在以下数据中，我使用了 quadprog 来求解线性模型。请注意，X1 只是截距。

X1 <- 1
X2 <- c(4.374, 2.3708, -7.3033, 12.0803, -0.4098, 53.0631, 13.1304, 7.3617, 16.6252, 27.3394)
X3 <- c(2.6423, 2.6284, 36.9398, 15.9278, 18.3124, 54.5039, 3.764, 19.0552, 25.4906, 13.0112)
X4 <- c(4.381, 3.144, 9.506, 15.329, 21.008, 38.091, 22.399, 13.223, 17.419, 19.87)
X <- as.matrix(cbind(X1,X2,X3,X4))
Y <- as.matrix(c(37.7,27.48,24.08,25.97,16.65,73.77,45.10,53.35,61.95,71.15))
M1 <- solve.QP(t(X) %*% X, t(Y) %*% X, matrix(0, 4, 0), c())$solution

挑战在于使某些系数受到约束。我知道我应该改变 Amet 和 bvac 参数（根据带有约束系数的线性回归）。但是，我不确定如何设置它，以便满足以下约束。

输出为 [1] 37.3468790 1.2872473 -0.0177749 -0.5988443，其中的值将是 X1、X2、X3 和 X4 的预测拟合值。

约束（受制于）……</p>

X2 <= .899

0 <= X3 <= .500

0 <= X4 <= .334

Answer 1

只是为了澄清一下：您给出了一个预期的输出“其中的值将被预测为 X1、X2、X3 和 X4 的拟合值”。我假设您指的是估计的（未拟合的）参数值。

在rstan. rstan是 的 R 接口Stan，而后者又是用于统计贝叶斯推理的概率编程语言。

这是基于您提供的示例数据的示例。

1. 首先，让我们定义我们的模型。请注意，我们限制参数beta2，beta3并beta4根据您的要求。

   model <- "
   data {
       int<lower=1> n;   // number of observations
       int<lower=0> k;   // number of parameters
       matrix[n, k] X;   // data
       vector[n] y;      // response
   }
   
   parameters {
       real beta1;                                                  // X1 unconstrained
       real<lower=negative_infinity(), upper=0.899> beta2;          // X2 <= .899
       real<lower=0, upper=0.5> beta3;                              // 0 <= X3 <= 0.5
       real<lower=0, upper=0.334> beta4;                            // 0 <= X4 <= 0.334
       real sigma;                                                  // residuals
   }
   
   model {
       // Likelihood
       y ~ normal(beta1 * X[, 1] + beta2 * X[, 2] + beta3 * X[, 3] + beta4 * X[, 4], sigma);
   }"
   

2. 接下来，我们将模型拟合到您的样本数据。

   library(rstan)
   rstan_options(auto_write = TRUE)
   options(mc.cores = parallel::detectCores())
   df <- cbind.data.frame(Y, X)
   fit <- stan(
       model_code = model,
       data = list(
           n = nrow(df),
           k = ncol(df[, -1]),
           X = df[, -1],
           y = df[, 1]))
   

3. 最后，让我们以整洁的方式打印参数估计

   library(broom)
   tidy(fit, conf.int = TRUE)
   ## A tibble: 5 x 5
   #  term  estimate std.error conf.low conf.high
   #  <chr>    <dbl>     <dbl>    <dbl>     <dbl>
   #1 beta1   29.2      6.53   16.9        42.8
   #2 beta2    0.609    0.234   0.0149      0.889
   #3 beta3    0.207    0.138   0.00909     0.479
   #4 beta4    0.164    0.0954  0.00780     0.326
   #5 sigma   16.2      5.16    9.42       29.1
   

   我们还可以绘制参数估计值，包括 CI。

请注意，参数估计与施加的约束一致。