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; iRows; 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.