R中具有等式和不等式约束系数的回归答案

【问题标题】：Regression with equality and inequality constrained coefficients in RR中具有等式和不等式约束系数的回归
【发布时间】：2018-07-24 09:32:42
【问题描述】：

我正在尝试使用 RSS 获得估计的约束系数。 beta 系数被限制在 [0,1] 和总和为 1 之间。此外，我的第三个参数被限制在 (-1,1) 之间。利用下面的方法，我可以使用模拟变量获得一个很好的解决方案，但是在我的真实数据集上实施该方法时，我一直在得出一个非唯一的解决方案。反过来，我想知道是否有一种在数值上更稳定的方法来获得我的估计参数。

set.seed(234)
k = 2
a = diff(c(0, sort(runif(k-1)), 1))
n = 1e4
x = matrix(rnorm(k*n), nc = k)
a2 = -0.5
y = a2 * (x %*% a) + rnorm(n)
f = function(u){sum((y - u[3] * (x %*% u[1:2]))^2)}
g = function(v){

v1 = v[1]
v2 = v[2]
u = vector(mode = "double", length = 3)

# ensure in (0,1)
v1 = 1 / (1 + exp(-v1))

# ensure add up to 1
u[1:2] = c(v1, 1 - sum(v1))

# ensure between [-1,1]
u[3] = (v2^2 - 1) / (v2^2 + 1)
u
}

res = optim(rnorm(2), function(v) f(g(v)), hessian = TRUE, method = "BFGS")
eigen(res$hessian)$values
res$convergence
rbind(Est = res$par, SE = sqrt(diag(solve(res$hessian))))
rbind(g(res$par),c(a,a2))

向http://zoonek.free.fr/blosxom/R/2012-06-01_Optimization.html致敬

【问题讨论】：

当stackoverflow.com/users/1129973/vincent-zoonekynd 说明htat 页面上的sum-to-one 约束时，他使用了quadprog。也许他会出现并帮助你。
@42- 很有帮助
@user5087936 请看下面我的示例，了解如何在 Stan/RStan 中实现具有约束参数的模型。了解模型在您的真实数据上的表现会很有趣。

标签： r optimization constraints regression

【解决方案1】：

由于到目前为止还没有直接回答您的问题，因此我想展示一种如何在 Stan/RStan 中实现参数约束模型的方法。您应该尝试使用您的真实数据。

进行贝叶斯推理的优势在于为您的（受约束的）模型参数提供后验概率。然后可以轻松计算包括置信区间在内的点估计值。

首先，我们加载库并将 RStan 设置为存储编译后的模型并使用多个内核（如果可用）。
```
library(rstan);
rstan_options(auto_write = TRUE);
options(mc.cores = parallel::detectCores());
```

我们现在定义我们的 Stan 模型。在这种情况下，它非常简单，我们可以将 RStan 的 simplex 数据类型用于总和为 1 的非负值向量。

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 a2;          // a2 is a free scaling parameter
    simplex[k] a;     // a is constrained to sum to 1
    real sigma;       // residuals
}

model {
    // Likelihood
    y ~ normal(a2 * (X * a), sigma);
}"

Stan 支持各种受约束的数据类型；我建议您在Stan manual 上大量了解更复杂的示例。

使用您原始问题中的样本数据，我们可以运行我们的模型：

# Sample data
set.seed(234);
k = 2;
a = diff(c(0, sort(runif(k-1)), 1));
n = 1e4;
x = matrix(rnorm(k * n), nc = k);
a2 = -0.5;
y = a2 * (x %*% a) + rnorm(n);

# Fit stan model
fit <- stan(
    model_code = model,
    data = list(
        n = n,
        k = k,
        X = x,
        y = as.numeric(y)),
    iter = 4000,
    chains = 4);

运行模型仅需几秒钟（在解析器内部翻译并用 C++ 编译模型之后），完整结果（所有参数的后验分布以数据为条件）存储在fit。

我们可以使用summary检查fit的内容：

# Extract parameter estimates
pars <- summary(fit)$summary;
pars;
#               mean      se_mean          sd          2.5%           25%
#a2       -0.4915289 1.970327e-04 0.014363398    -0.5194985    -0.5011471
#a[1]      0.7640606 2.273282e-04 0.016348488     0.7327691     0.7527457
#a[2]      0.2359394 2.273282e-04 0.016348488     0.2040952     0.2248482
#sigma     1.0048695 8.746869e-05 0.007048116     0.9909698     1.0001889
#lp__  -5048.4273105 1.881305e-02 1.204892294 -5051.4871931 -5048.9800451
#                50%           75%         97.5%    n_eff      Rhat
#a2       -0.4916061    -0.4819086    -0.4625947 5314.196 1.0000947
#a[1]      0.7638723     0.7751518     0.7959048 5171.881 0.9997468
#a[2]      0.2361277     0.2472543     0.2672309 5171.881 0.9997468
#sigma     1.0048994     1.0095420     1.0187554 6492.930 0.9998086
#lp__  -5048.1238783 -5047.5409682 -5047.0355381 4101.832 1.0012841

你可以看到a[1]+a[2]=1。

绘制包括置信区间的参数估计也很容易：

plot(fit);

【讨论】：

【解决方案2】：

解决具有等式和不等式约束的优化问题的最简单方法很可能是通过“增广拉格朗日”方法。例如，在 R 中，这是在 alabama 包中实现的。

# function and gradient
fn = function(u){sum((y - u[3] * (x %*% u[1:2]))^2)}
gr = function(u) numDeriv::grad(fn, u)

# constraint sum(u) == 1
heq = function(u) sum(u) - 1
# constraints 0 <= u[1],u[2] <= 1; -1 <= u[3] <= 1
hin = function(u) c(u[1], u[2], 1-u[1], 1-u[2], u[3]+1, 1-u[3])

sol_a = alabama::auglag(c(0.5, 0.5, 0), fn, gr, hin=hin, heq=heq)
sol_a
## $par
## [1]  1.0000000  0.3642904 -0.3642904
## $value
## [1] 10094.74
## ...
## $hessian
##             [,1]        [,2]        [,3]
## [1,] 15009565054  9999999977  9999992926
## [2,]  9999999977 10000002578  9999997167
## [3,]  9999992926  9999997167 10000022569

对于包含“增强拉格朗日”过程的其他包，请参阅 CRAN 任务视图优化。

【讨论】：