Skip to main content

Case Studies

Case study: Thermal Expansion of Copper

nistlogo 
NIST info: This case study illustrates the use of a class of nonlinear models called rational function models. The data set used is the thermal expansion of copper related to temperature. This data set was provided by the NIST scientist Thomas Hahn. Full description can be found here.

Data: The response variable for this data set is the coefficient of thermal expansion for copper. The predictor variable is temperature in degrees kelvin. There were 236 data points collected (data avaiable here).

To perform non-linear regression on collected data, we willl use Dew Stats tMtxNonLinReg control. For non-linear fit, we have to define the fit function, initial values for fit parameters and (optionaly) optimization method used in nonlinear regression. We begin with defining the quadratic-quadratic rational function:

project d32 
function FracQuadQuad(const B: Array of TSample; X: TSample): TSample; 
begin 
   Result := (B[0] + X*B[1] + X*X*B[2])/(1.0 + B[3]*X + B[4]*X*X);
end;
project cpp 
double FracQuadQuad(const double* B, const int param_size, double X);
{
   return (B[0]+X*B[1]+X*X*B[2])/(1.0 + B[3]*X + B[4]*X*X);
}
project cs 
double FracQuadQuad(double[] B, double X)
{
 return (B[0]+X*B[1]+X*X*B[2])/(1.0 + B[3]*X + B[4]*X*X);
}


This is the function we'll be fitting to measured data. Next, we need initial estimates for the regression coefficients i.e. the values for B[0]..B[4]. As suggested at NIST pages, we'll use linear fit on selected points to get initial estimates for regression coefficients. For this we can use the MulLinRegress routine. Here is the implementation:

project d32 
// Do an initial estimate on fractional model by performing a
// linear fit on y = B[0] + B[1]*x + ... + B[nom]*x^nom
// - B[nom+1]*x*y - B[nom]*x*x*y - ... - B[nom+denom]*x^denom
procedure FracInitEstimates(x,y: TVec; b: TVec; nom,denom: Integer);
   var m1,m2,m3: Matrix;
   i: Integer;
begin
   if (nom+denom+1)<>(x.Length) then Raise Exception.Create('(nom+denom+1)<>(x.Length)');
   // ..x*x, x, 1
   m1.VanderMonde(nom+1,x);
   // flip
   m1.FlipHor;
   // ..x*x, x, 1
   m2.VanderMonde(denom+1,x);
   // flip and remove 1
   m2.Resize(m2.Rows,denom);
   m2.FlipHor;
   // multiply each row with y[i]
   for i := 0 to m2.Rows - 1 do m2.Mul(y.Values[i],i*m2.Cols,m2.Cols);
   m3.ConcatHorz([m1,m2]);
   MulLinRegress(y,m3,b,false);
end;
project cpp 

void FracInitEstimates(Vector x, Vector y, Vector b, int nom, int denom);
{
   Matrix m1,m2,m3;
   try
   m1->VanderMonde(nom+1,x);
   m1->FlipHor();
   m2->VanderMonde(denom+1,x);
   m2->Resize(m2->Rows,denom);
   m2->FlipHor();
   // multiply each row with y[i]
   for (int i=0; iRows; i++) m2(i) *= y[i];
 
   TMtx* conc[2];
   conc[0] = m1;
   conc[1] = m2;
   m3->ConcatHorz(conc,1);
   MulLinRegress(y,m3,b,false);
}
project cs 
using Dew.Stats.Units;
using Dew.Stats;
using Dew.Math.Units;
using Dew.Math;

namespace Dew.Tests
{
   // Do an initial estimate on fractional model by performing a
   // linear fit on y = B[0] + B[1]*x + ... + B[nom]*x^nom
   // - B[nom+1]*x*y - B[nom]*x*x*y - ... - B[nom+denom]*x^denom
   private void FracInitEstimates(TVec x, TVec y, TVec b, int nom, int denom)
   {
      Matrix m1=new Matrix(0,0);
      Matrix m2=new Matrix(0,0);
      Matrix m3=new Matrix(0,0);
      // ..., x*x, x, 1
      m1.VanderMonde(nom+1,x);
      // flip
      m1.FlipHor();
      // ..., x*x, x, 1
      m2.VanderMonde(denom+1,x);
      // flip and remove 1
      m2.Resize(m2.Rows,denom);
      m2.FlipHor();
      // multiply each row with y[i]
      for (int i=0; i<m2.Rows; i++)
         m2.Mul(y.Values[i],i*m2.Cols,m2.Cols);
      m3.ConcatHorz(m1,m2);
      Regress.MulLinRegress(y,m3,b,false,null,null);
   }
}

Now all we must do is put all this together and perform the non-linear fit. The call to FitModel() routine first loads data, then performs initial estimation of regression coefficients and then uses initial estimates in actual non-linear regression:

project d32 
procedure FitModel;
var nlr: TMtxNonLinReg; 
   Residuals, initx, inity: Vector;
   ResSTdDev: Tsample;
begin
   nlr := TMtxNonLinReg.Create(nil);
   try
      // load data
      nlr.X.LoadFromFile('nist_001_temp.vec');
      nlr.Y.LoadFromFile('nist_001_coeff.vec');
      // Define and connect to regression function
      nlr.RegressFunction := FracQuadQuad;
 
      // initial values for coefficients
      initx.SetIt(false,[10,50,120,200,800]);
      inity.SetIt(false,[0,5,12,15,20]);
      FracInitEstimates(initx,inity,nlr.B,2,2);
      // nlr.B returns 0.301e+1,0.369,-0.683e-2,-0.112e-1,-0.306e-3]
 
      // Perform non-linear regression
      nlr.Recalc;
      // Calculate residuals
      Residuals.Sub(nlr.YCalc,nlr.Y);
   finally
      nlr.Free;
   end;
end;>
project cpp 
void FitModel(void);
{
   Vector residuals, initx, inity;
   TMtxNonLinReg* nlr = new TMtxNonLinReg(NULL);
   try
   {
      nlr->LoadFromFile("nist_001_temp.vec");
      nlr->LoadFromFile("nist_001_coeff.vec");
      // connect to regression function
      nlr->RegresFun = FracQuadQuad;

      // initial estimates for coefficients
      initx->SetIt(OPENARRAY(TSample,(10,50,120,200,800)));
      inity->SetIt(OPENARRAY(TSample,(0,5,12,15,20)));
      FracInitEstimates(initx,inity,nlr->B,2,2);

      // perform non-linear regression
      nlr->Recalc();
      residuals->Sub(nlr->YCalc,nlr->Y);
   }
   __finally
   {
      nlr->Free();
   }
}
project cs 
namespace Dew.Tests
{
   private void FitModel()
   {
      TMtxNonLinReg nlr = new TMtxNonLinReg(this);
      Vector residuals = new Vector(0);
      Vector initx = new TVector(0);
      Vector inity = new Vector(0);
      try
      {
         // load data
         nlr.X.LoadFromFile("nist_001_temp.vec");
         nlr.Y.LoadFromFile("nist_001_coeff.vec");
         // connect to regression function
         nlr.RegressFun += new Dew.Stats.TRegressFun(FracQuadQuad);

         // initial values for coefficients
         initx.SetIt(new double[] {10,50,120,200,800});
         inity.SetIt(new double[] {0,5,12,15,20});
         FracInitEstimates(initx,inity,nlr.B,2,2);
         // nlr.B returns 0.301e+1,0.369,-0.683e-2,-0.112e-1,-0.306e-3]
 
         // Perform non-linear regression
         nlr.Recalc();
         // Calculate residuals
         residuals.Sub(nlr.YCalc,nlr.Y);
      }
      finally
      {
         nlr.Free();
      }
   }
}

Using these results we can then calculate residuals standard deviation, residuals sum of squares and additional statistics.

Additional things to try :

  • Change optimization method (tMtxNonLinReg.OptimizationMethod property).
  • Instead of quadratic rational function, try fitting cubic rational function. Check NIST page which temperature and coefficient values to use for initial estimates.
  • Try using the RegModels.FracFit routine to perform fit).
  • Plot data, fitted data and residuals by using DrawValues routine.
  • Compare the results with quoted results at NIST pages.