Curve Fitting: Fitting a Curve to a Set of Data Points

The following PROC OPTMODEL statements declare an index set and parameters and then read the input data:

proc optmodel;
   set POINTS;
   num x {POINTS};
   num y {POINTS};
   read data xy_data into POINTS=[_N_] x y;

The following NUM statement declares the order parameter, which is later set to 1 for Problems (1) and (2) and set to 2 for Problem (3):

   num order;
   var Beta {0..order};
   impvar Estimate {i in POINTS}
      = Beta[0] + sum {k in 1..order} Beta[k] * x[i]^k;

The following statements encode the linearization of the norm:

   var Surplus {POINTS} >= 0;
   var Slack {POINTS} >= 0;
   min Objective1 = sum {i in POINTS} (Surplus[i] + Slack[i]);
   con Abs_dev_con {i in POINTS}:
      Estimate[i] - Surplus[i] + Slack[i] = y[i];

The following statements (which are not used but are shown here for comparison) encode an alternative linearization of the norm that requires half as many variables but twice as many constraints:

   var AbsDeviation {POINTS} >= 0;
   min Objective1 = sum {i in POINTS} AbsDeviation[i];
   con Abs_dev_con1 {i in POINTS}:
      AbsDeviation[i] >= Estimate[i] - y[i];
   con Abs_dev_con2 {i in POINTS}:
      AbsDeviation[i] >= y[i] - Estimate[i];

The following additional declarations encode the linearization of the norm:

   var MinMax;
   min Objective2 = MinMax;
   con MinMax_con {i in POINTS}:
      MinMax >= Surplus[i] + Slack[i];

The following statements (which are not used but match the formulation in Williams (1999)) encode an alternative linearization of the norm that requires the same number of variables but twice as many constraints:

   var MinMax;
   min Objective2 = MinMax;
   con MinMax_con1 {i in POINTS}:
      MinMax >= Surplus[i];
   con MinMax_con2 {i in POINTS}:
      MinMax >= Slack[i];

The following NUM statements use the ABS function, the .sol variable suffix, and the MAX aggregation operator to compute the two norms from the optimal solution:

   num sum_abs_dev = sum {i in POINTS} abs(Estimate[i].sol - y[i]);
   num max_abs_dev = max {i in POINTS} abs(Estimate[i].sol - y[i]);

The following PROBLEM statement specifies the variables, objective, and constraints for minimization:

   problem L1 include
      Beta Surplus Slack
      Objective1
      Abs_dev_con;

The following PROBLEM statement specifies the additional variables, objective, and constraints for minimization:

   problem Linf from L1 include
      MinMax
      Objective2
      MinMax_con;

The following statements specify a straight line fit, switch the focus to problem L1, solve Problem (1), print the results, and store the y-values that are predicted by the optimal solution:

   order = 1;
   use problem L1;
   solve;
   print sum_abs_dev max_abs_dev;
   print Beta;
   print x y Estimate Surplus Slack;
   create data sol_data1 from [POINTS] x y Estimate;

Figure 11.1 shows the output from the linear programming solver for Problem (1).

The following statements switch the focus to problem Linf, solve Problem (2), print the results, and store the y-values that are predicted by the optimal solution:

   use problem Linf;
   solve;
   print sum_abs_dev max_abs_dev;
   print Beta;
   print x y Estimate Surplus Slack;
   create data sol_data2 from [POINTS] x y Estimate;

Figure 11.2 shows the output from the linear programming solver for Problem (2).

The following statements specify a quadratic curve fit, solve both parts of Problem (3), print the results, and store the y-values that are predicted by the optimal solutions:

   order = 2;
   use problem L1;
   solve;
   print sum_abs_dev max_abs_dev;
   print Beta;
   print x y Estimate Surplus Slack;
   create data sol_data3 from [POINTS] x y Estimate;

   use problem Linf;
   solve;
   print sum_abs_dev max_abs_dev;
   print Beta;
   print x y Estimate Surplus Slack;
   create data sol_data4 from [POINTS] x y Estimate;
quit;

Figure 11.3 shows the output from the linear programming solver for the first part of Problem (3).

Figure 11.4 shows the output from the linear programming solver for the second part of Problem (3).

You can find a higher-order polynomial fit simply by increasing the value of the order parameter. The dimensions of the Beta and Estimate variables are automatically updated when order changes.

The following PROC SGPLOT statements use the output data sets that are created by PROC OPTMODEL to display the results from Problems (1) and (2) in one plot:

data plot1;
   merge sol_data1(rename=(Estimate=Line1)) sol_data2(rename=(Estimate=Line2));
run;

proc sgplot data=plot1;
   scatter x=x y=y;
   series x=x y=Line1 / curvelabel;
   series x=x y=Line2 / curvelabel;
run;

Figure 11.5 shows the regression lines for Problems (1) and (2), as on page 319 of Williams (1999).

Figure 11.5: Regression Lines for Problems (1) and (2)

The following PROC SGPLOT statements use the output data sets that are created by PROC OPTMODEL to display the results from both parts of Problem (3) in one plot:

data plot2;
   merge sol_data3(rename=(Estimate=Curve1))
      sol_data4(rename=(Estimate=Curve2));
run;

proc sgplot data=plot2;
   scatter x=x y=y;
   series x=x y=Curve1 / curvelabel;
   series x=x y=Curve2 / curvelabel;
run;

Figure 11.6 shows the regression curves for Problem (3), as on page 319 of Williams (1999).

Figure 11.6: Regression Curves for Problem (3)