Dew Research - NIST case 01

Case study: Thermal Expansion of Copper

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:

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;

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);
}

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:

// 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;

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; i<m2->Rows; i++)  m2(i) *= y[i];
    
  TMtx* conc[2];
  conc[0] = m1;
  conc[1] = m2;
  m3->ConcatHorz(conc,1);
  MulLinRegress(y,m3,b,false);
}

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:

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;

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();
  }
}

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.

Back to Cases

Delphi sources

C++ Builder sources

C# sources