PROC OPTMODEL help needed - optimising a linear equation with restrictions

[Jonathan_A5]

I'd really like to be able to solve the following problem:

Given observed data x(1)....x(n) and a known fixed 'target' T, solve for parameters a0, a1, & a2, which minimise:

{ T - sum[i=1 to n] exp(a0+a1*x(i)+a2*x(i)^2)*x(i)  }^2

with the restriction that the sum of the exp(a0+a1*x(i)+a2*x(i)^2) terms equals n,

i.e. the mean of the exp(a0+a1*x(i)+a2*x(i)^2) terms equals 1.

I've tried PROC MODEL and PROC NLIN in SAS v9.4 but with no luck.

I'm now attempting to try with PROC OPTMODEL

My OPTMODEL test code is as follows:

* Create data with known parameters;
data d1;
do i = 1 to 100;
x=ranuni(4245);
w1=exp(1.0253056-3*x+0.4*x**2);* Weighting variable;
output;
end;
run;
proc means data=d1 n mean;var w1;run; * Mean weight of 1;
proc means data=d1 n mean;var x;run;
proc means data=d1 n mean;var x;weight w1;run; * Target (T)=0.2736;

proc optmodel;
set OBS;
num x{OBS};
read data d1 into OBS = [i] x ;
var a0 init 1, a1 init -2, a2 init 0.5;
con restrict: sum {i in OBS} exp(a0 + a1*x[i] + a2*x[i]**2) = sum{i in OBS} 1;
min objective = (0.2736 - sum {i in OBS} (exp(a0 + a1*x[i] + a2*x[i]**2)*x[i]))**2;
solve with nlp / tech=ip ms seed=21 MSMAXSTARTS=5000 msbndrange=50;
print a0 a1 a2;
quit;run;

The model does not converge anywhere near the a0, a1 & a2 values used in the 'd1' data step, even when initiated at values very close to those actually used.

What am I doing wrong? I suspect it's the way the objective function is specified....but just a hunch.

Any help gratefully received!

[RobPratt]

Your current objective is trying to force the mean and the sum to be equal.  You instead want to replace 0.2736 with 0.2736*card(OBS).

Here is an alternative approach that avoids hard-coding the target and also avoids duplication of the exponential expression:

proc optmodel;
   set OBS;
   num x{OBS};
   num w1{OBS};
   read data d1 into OBS = [i] x w1;
   num target = sum {i in OBS} w1[i]*x[i];
   num aInit {0..2} = [1, -2, 0.5];
   var a {j in 0..2} init aInit[j];
   impvar w {i in OBS} = exp(sum {j in 0..2} a[j]*x[i]^j);
   con restrict: sum {i in OBS} w[i] = card(OBS);
   min objective = (target - sum {i in OBS} w[i]*x[i])**2;
   solve with nlp / tech=ip ms seed=21 MSMAXSTARTS=5000 msbndrange=50;
   print aInit a;
quit;

[Jonathan_A5]

Hi Rob,

Many thanks for your suggestion. 

Please note that the w(i) in the problem are actually unknown - I just simulated them to see of the the model could converge on the three parameters (a1, s2, a3) that are used to calculate the w(i) in the dataset d1.

The only knowns in a 'real-life' situation will be the x(i) and the target.

The only restrictions are that the w(i) need to be positive (hence the use of the exponential function) and that the mean of the w(i) needs to be 1. 

I'm guessing the above would impact upon your suggested code. Happy to hear your thoughts.

Kind regards,

Jonathan

[RobPratt · Accepted solution]

Glad to help.  The only place I used the actual values of w1(i) was to compute the target as the weighted sum.  You could pass that single number in as an input instead:

proc sql;
   select sum(w1*x) into :target
   from d1;
quit;

proc optmodel;
   set OBS;
   num x{OBS};
   read data d1 into OBS = [i] x;
   num aInit {0..2} = [1, -2, 0.5];
   var a {j in 0..2} init aInit[j];
   impvar w {i in OBS} = exp(sum {j in 0..2} a[j]*x[i]^j);
   con restrict: sum {i in OBS} w[i] = card(OBS);
   min objective = (&target - sum {i in OBS} w[i]*x[i])**2;
   solve with nlp / tech=ip ms seed=21 MSMAXSTARTS=5000 msbndrange=50;
   print aInit a;
   print x w;
quit;

[Jonathan_A5]

Many thanks Rob!!

This was a great help. I've been a SAS user for 25+ years, but I've only just started using Proc OPTMODEL and there is the usual, rather steep, learning curve. It never ceases to amaze me the breadth and depth of modules within SAS...the more you know, the less you know, etc.

best wishes and kind regards,

Jonathan