使用 Quadprog 对系数 Ranges/Limits 的线性模型的约束

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 的预测拟合值。

约束（受制于）…

X2 <= .899

0 <= X3 <= .500

0 <= X4 <= .334

Answer 1

澄清一下：您给出了预期的输出 "where the values would be predicted fitted values of X1, X2, X3 and X4"。我假设您指的是 estimated（未拟合）参数值。

在 rstan 中对数据建模时，实现约束参数非常简单。 rstan 是 Stan 的 R 接口，后者又是一种用于统计贝叶斯推理的概率编程语言。

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

首先，让我们定义我们的模型。请注意，我们根据您的要求限制参数 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);
}"

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

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]))

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

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。

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

使用 Quadprog 对系数 Ranges/Limits 的线性模型的约束

Constraints for a Linear Model Using Quadprog for Ranges/Limits on Coefficients

r

linear-regression

quadprog