How to implement positivity constraints to a function of data and parameters

Modeling
1

Dear Stan comunity,

I’m still new to Stan, and I’ve been facing some problems while writing my code.
I am trying to estimate four parameters, namely C_1, C_2, C_3 and C_4. They are the coefficients of 3rd degree function f(x) = C_1 x^3 + C_2 x^2 + C_3 x + C_4. The only prior information I have about the coefficients is C_4 > 0. Additionally, I know that f(x) > 0 in the range d_{\text{min}} < x < d_{\text{max}}. These values d_{\text{min}} and d_{\text{max}} represent, respectively, the minimum and maximum values of my data. In my problem, x is always between d_{\text{min}} and d_{\text{max}}, so I don’t have to worry about what happens outside this interval.

So far my code is similar to this:

functions {

  real f(real x, real C_1, real C_2, real C_3, real C_4) {
  real my_f = C_1 * pow(x, 3) + C_2 * pow(x, 2) + C_3 * x + C_4;
  return my_f;
  }

  vector complicated_function(C_1, C_2, C_3, C_4) {
  real x;
  // Do some calculation to find "x", which is a function of C_1, C_2, C_3 and C_4
  real y = f(x, C_1, C_2, C_3, C_4);
  real z;
  // Do some calculation to find "z", which is a function of "y".
  return z;
  }
}

data {
  int n;
  vector[n] d;
}

transformed data {
  real d_min = min(d);
  real d_max = max(d);
}

parameters {
  real C1;
  real C2;
  real C3;
  real<lower=0> C4;
}

transformed parameters {
  vector[n] ff = complicated_function(C_1, C_2, C_3, C_4);
}

model {
  d ~ normal(ff, 1); // Just an example.
}

I got stuck when I tried to implement the constraint C_1 x^3 + C_2 x^2 + C_3 x + C_4 > 0 for d_{\text{min}} < x < d_{\text{max}}. I have the impression that I should combine the blocks “transformed data” and “transformed parameters”, but I am not sure how. I would really appreciate if someone could help me with this.

Thanks in advance.

Best regards,
Rodrigo

2

Hi Rodrigo,

I have done some work in a paper that could be applied here.

Without getting into too much technical details I essentially penalize the observable outcome (in your case the sum of the coefficients). From your question you want f(x) or ff (in your code) >0 (I’m not really sure where x comes from in your example but for my explanation it doesn’t matter).

target += uniform_lpdf(ff |0, 10000) #from 0 to some very large number

Now Stan will reject any values of ff less than 0 but you will have to specify valid initial values otherwise the loglikelihood will be negative infinity.

I’d be interested in understanding the context of your regression i.e. what biological/economic model requires the function to be greater than zero? It would be great to add some further potential use cases to the paper I linked (that is currently undergoing peer review) :)

Philip

2 Likes
3

Hi Philip,

Thank you very much for your reply. I managed to solve my problem. My solution was a bit similar to the strategy you recommended. I’ll omit the solution here just to avoid a lot of technical details that might not be relevant for you, but let me know in case you’re interested to know how I solved it.

My problem is actually related to heat transfer problems. Based on temperature measurements, I estimated the thermal conductivity of a material. In my problem, the thermal conductivity k is a function of the temperature T, so k = k(T). This property is strictly a positive number, and therefore I need to make sure that k(T) > 0.

Best,
Rodrigo

4

Thanks @Rodrigo_Silva, makes sense.
Glad you got your problem solved.

Philip

1 Like