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

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.