具有相等和不等式的回归约束了R中的系数
Posted
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了具有相等和不等式的回归约束了R中的系数相关的知识,希望对你有一定的参考价值。
我试图使用RSS获得估计的约束系数。 β系数被约束在[0,1]和sum之间。另外,我的第三个参数被约束在(-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致敬
答案
由于到目前为止还没有直接回答你的问题,我想说明如何在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
中。 - 我们可以使用
fit
检查summary
的内容:# 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);
另一答案
解决具有相等和不等式约束的优化问题的最简单方法很可能是通过“增广拉格朗日”方法。在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任务视图。
以上是关于具有相等和不等式的回归约束了R中的系数的主要内容,如果未能解决你的问题,请参考以下文章