Inequality constraint aX+b>0 in the sigma of normal distribution

Modeling
1

Hi, everyone. I would like to use inequality constraint like Y \sim normal(mu, aX+b), where aX+b > 0 and X is data. However, in most samplings with different tried codes, I had the errors “Scale parameter is negative value, but must be > 0!” and finally “Stan model ‘model name’ does not contain samples.” The basic structure of my Stan code is as follows:

data {
  int N; // number of samples
  int Group[N];
  int NGroup; // number of groups
  real Y[N];
  real X[N]; // X > 0
}

parameters {
  real mu[NGroup];
  real sa[NGroup];
  real sb[NGroup];
}

model {
  sa ~ student_t(4, 0, 100);
  sb ~ student_t(4, 0, 100);
  for (n = 1 : N) {
    if (sa[Group[n]] * X[n] + sb[Group[n]] < 0)
      target += negative_infinity();
    Y[n] ~ normal(mu[Group[n]], sa[Group[n]] + sb[Group[n]] * X[n]); 
  }
}

I had the “not contain samples” error even when I wrote the “reject” statement instead of “negative_infinity.” Since aX+b > 0 means b > -aX, I also tried setting in the parameters block:

  real sa[NGroup];
  real<lower=-sa[NGroup]*X[N]> sb[NGroup];

However, I still got the initially mentioned errors. The same applies to the case where the transformed parameters block has the constraint:

transformed parameters {
  real<lower=0> s[NGroup];
  for (n in 1 : N)
    s[Group[n]] = sa[Group[n]] * X[n] + sb[Group[n]];
}

If I set these constraints both in the parameters and transformed parameters, sampling ran without errors, but the posterior sa was almost zero (e.g. 0.01 in median), although X obviously affects the variance of Y in the measured data.

Could you suggest any solution to this? I perform sampling in R 3.5.1 with the “rstan” package; Stan version is 2.18.1. Thank you in advance for your kind help.

2

The easiest thing to do is make sigma be some positive transformation of a real number like exp(a * X + b). Perhaps the second easiest is to declare sigma as a positive parameter and parameterize its conditional mean.

The hardest thing to do is constrain b so that aX + b is positive for all values of X. You can technically do that with something like

data {
  int N; // number of samples
  int Group[N];
  int NGroup; // number of groups
  vector[N] Y;
  vector<lower=0>[N] X; 
}

parameters {
  vector[NGroup] mu;
  vector[[NGroup] sa;
  vector<lower = -min(sa[Group] .* X)>[NGroup] sb;
}
model {
  sa ~ student_t(4, 0, ???);  // don't use 100 for the scale
  sb ~ ??? // don't use half-t
  Y ~ normal(mu[Group], sa[Group] + sb[Group] .* X); // don't use a loop
}

but it is hard to come up with a reasonable prior on sb over the positive real numbers in light of sa.

3

Thank you so much, Ben. I set exp(aX+b) for sigma and succeeded in sampling. Your other comments about priors and syntax are also very helpful.