每当我有一个复杂的模型要拟合时，我通常只是直接拟合它，rstan因为它非常适合拟合高度约束的系数，并且因为它很容易包含变量的惩罚和转换。即使我没有明确拟合贝叶斯模型也是如此。

这就是我为您的特定问题所做的工作。

library(rstan)

set.seed(1880)

N       <- 1500
d       <- c(1/2, 2/pi, 2/3)
x       <- c(2, 1, 3)
limit   <- 5
d%*%x <= limit
> TRUE
A       <- cbind(1, rnorm(N), rnorm(N))
b.hat   <- A%*%x
tau     <- 5
wgt     <- rexp(N)
Sigma   <- tau*wgt
b       <- rnorm(N, mean=b.hat, sd=Sigma)

constrained.reg <- "
    data{
        int<lower=1>        N;
        int<lower=1>        K;
        vector<lower=0>[N]  wgt;
        matrix[N,K]         A;
        vector[N]       b;
        vector[K]       d;
        real            limit; // s.t. d*x<=limit
    }
parameters{
    real<upper=limit>   c; // this is the largest possible value of x%*%d.
    simplex[K]      sim_x;
    real<lower=0>       tau;
}
transformed parameters {
    vector[K]   x;
    vector[N]   b_hat;
    vector[N]   Sigma;

    x       <- d .*sim_x /c;
    b_hat   <- A*x;
    Sigma   <- tau*wgt;
}
    model{
        b ~ normal(b_hat, Sigma);
        increment_log_prob(-2*log(tau)); // uniform prior on beta, noninformative prior on tau
    }
    generated quantities{
        vector[N]   resid;
        resid   <- (b_hat-b) ./Sigma;
    }
"
fake.data   <- list(N=N, A=A, K=3, b=b, wgt=wgt, d=d, limit=limit)

fit.test    <- stan(model_code=constrained.reg, data=fake.data, iter=10)

system.time(fit     <- stan(fit=fit.test, iter=1000, data=fake.data))
print(fit, c("x", "tau")); x

我意识到我很密集，我们可以通过采样一个与最大允许点积结果一样大的值然后进行适当的转换来强制执行不等式。

     mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
x[1] 1.99       0 0.01 1.98 1.98 1.99 1.99  2.00  1645 1.00
x[2] 0.99       0 0.01 0.97 0.98 0.99 0.99  1.00   624 1.00
x[3] 3.00       0 0.01 2.98 2.99 3.00 3.01  3.02   945 1.00
tau  4.82       0 0.09 4.62 4.76 4.82 4.88  5.00   558 1.01

这些结果在我看来很好。

您正在寻找 mgcv 包。使用我们之前使用的玩具数据，它工作得很好。（我不确定为什么rstan对它的结果如此有信心......我还在研究它。）

set.seed(1880)

N       <- 1500
d       <- c(1/2, 2/pi, 2/3)
x       <- c(2, 1, 3)
limit   <- 5
d%*%x <= limit

A       <- cbind(1, rnorm(N), rnorm(N))
b.hat   <- A%*%x
wgt     <- rexp(N)
b       <- rnorm(N, mean=b.hat, sd=wgt)

library(mgcv)

pin <- c(1.5, .75, 2.5)
Ain <- matrix(d, nrow=1)

M   <- list(y=b, w=wgt, X=A, p=pin, Ain=-Ain, bin=-limit, C=matrix(1, ncol=0, nrow=0))
pcls(M)

1.8844996 0.9421333 2.9770852

默认情况下，这个包中的不等式被翻转到另一个方向。所以我们必须将两边都乘以。$-1$