LBFGS Optimization to maximize function with inequality constraints on sum

Modeling
, , ,
2

Bounds aren’t calculated in the transformed parameters block, they are only enforced. To get Stan to enforce the bounds on the total you need to calculate each bound separately given the previous. You don’t need a uniform distribution prior on the total. With the code below it will automatically enforce that. Since there’s not additional data though I believe this optimization will always have the total be equal to the max.

Also, Stan will complain if the bound is included so we need to add a small amount to allow that.

This works where the bounds are calculated for each rate given the estimate of the previous rate.

functions {   
   real lb_ub_lp (real y, real lb, real ub) {
    target += log(ub - lb) + log_inv_logit(y) + log1m_inv_logit(y);
    
    return lb + (ub - lb) * inv_logit(y);
  }
}
data{
  // shape
  int<lower=0> N;
  
  // data
  real<lower=0> mult;
  vector[N] x;
  
  // constraints
  real<lower=0> max_value;
}
transformed data {
  vector[N] x_mult = x * mult;
}
parameters{
  vector[N] rates_raw;
}

transformed parameters{
  vector[N] rates;
  vector[N] x_rate_mult;
  
  rates[1] = lb_ub_lp(rates_raw[1], 0, max_value / x_mult[1]);
  x_rate_mult[1] = rates[1] * x_mult[1];
  {
    real max_ = max_value;
    for (i in 2:N) {
      max_ += -rates[i - 1] * x_mult[i - 1];
      rates[i] = lb_ub_lp(rates_raw[i], 0, max_ / x_mult[i]);
      x_rate_mult[i] = rates[i] * x_mult[i];
    }
  }
  real<lower=0, upper=max_value + 1> total = dot_product(rates, x_mult);
}
model{
  target += x_rate_mult;
  
}
1 Like