ECE 3340

Part A - Model Fit and Regression

1. Linear regression introduction

Use the APP below to follow the lecture (copy and paste into a new notebook)

1.1 Illustration

Suppose we have a known relationship between two variables: 
                 y= a x+b 
Suppose we have data like this (you can also use the APP above):

RandomSeed[Date[]];
Manipulate[data=Table[{i, 2*i+5+RandomVariate[NormalDistribution[0.,σnoise]]},{i,1,20}] ;
datplot1P=Graphics[{Hue[0.,1,1.],PointSize[0.03],Point /@ data}];
datplot1=Show[ListLinePlot[data, PlotStyle-> Hue[0.,1.,1.], PlotRange-> All], datplot1P]
,{{σnoise,4.},0.01,10}]

How do we find coefficients a and b?

1.2 Linear regression

One way to do this is that we choose a and b such that (OverHat[(Subscript[y, i])]) = a (Subscript[x, i]) + b is as close to the data value (Subscript[y, i]) as possible. But we have many points, which one should be close? One way to decide how close our fit (OverHat[(Subscript[y, i])]) = a (Subscript[x, i]) + b  to the real values (Subscript[y, i]) is to make this measurement:
                          (∑_(i=1)^n)((OverHat[(Subscript[y, i])]) - (Subscript[y, i]))^2 as small as possible.
Minimizing the function:             (∑_(i=1)^n)(a (Subscript[x, i]) + b-(Subscript[y, i]))2  
leads to expression for a and b:
                  (OverHat[a])=     (((((∑_(i=1)^n)(Subscript[x, i])(Subscript[y, i]))/n)  -(((∑_(i=1)^n)(Subscript[x, i]) (∑_(i=1)^n)(Subscript[y, i]))/(n^(2  ))))/((((∑_(i=1)^n)((Subscript[x, i])^2))/n)-((((((∑_(i=1)^n)(Subscript[x, i]))/n)))^2)))=(((OverBar[x y])  - (OverBar[x])  (Overscript[y, _ ]))/((OverBar[(x^2)]) - ((OverBar[x])^2)))
                  (OverHat[b]) = (OverBar[y])-(OverHat[a])(OverBar[x])   

This is called linear regression and the method is known as “least square.” Notations (OverHat[a]) and (OverHat[b]) mean that they are only estimates of the true (but unknown) values of a and b. For each data series, we have different  (OverHat[a]) and (OverHat[b]). Thus, the estimates themselves are statistical variables.
(Note, we use the "bar" symbol to define the mean of a data series. For example,
                        (OverBar[x]) =((((Subscript[x, 1])+(Subscript[x, 2])+...+(Subscript[x, n])))/n)

One can use Mathematica function Fit or LinearModelFit

x=.;
model=Fit[ data,{1,x},{x}]

Show[Plot[model,{x,0,20},PlotStyle->{RGBColor[0.,0.,1]}],datplot1P]

If we use the formula, here is the way:

x=.;
data=Table[{i, 2*i+5+RandomVariate[NormalDistribution[0.,4.]]},{i,1,20}] ;
model=Fit[ data,{1,x},{x}]; {bF,aF}={model[[1]],model[[2,1]]};
{xdat1,ydat1}=Transpose[data] ;
a=((Mean[xdat1*ydat1] -Mean[xdat1]*Mean[ydat1])/(Mean[xdat1^2]-Mean[xdat1]^2));
b=Mean[ydat1] - a*Mean[xdat1];
Grid[Join[{{Style["Linear regression coef",Bold],"calc.","Fit"}}
,MapThread[Prepend,{{{a,aF},{b,bF}}
,{"a (slope)","b (intercept)"}}]]
,Frame-> All,BaseStyle->{20,FontFamily-> "Calibri"}
,Background->{{Hue[0.16,1,1],},{Hue[0.4,1,1],}}]

1.3 Correlation: how good is the “fit”?

Consider this example:

RandomSeed[Date[]];
Manipulate[noisydata=Table[{i, 2*i+5+RandomVariate[NormalDistribution[0.,σnoise]]},{i,1,20}] ;
{xdat2,ydat2}=Transpose[noisydata] ;
datplot2P=Graphics[{Hue[0.65,1,1.],PointSize[0.03],Point /@ noisydata}];
datplot2=Show[ListLinePlot[noisydata, PlotStyle-> Hue[0.6,1.,1.], PlotRange-> All], datplot2P];
Show[datplot2,datplot1]
,{{σnoise,15.},5,20}]

noisydata=Table[{i, 2*i+5+RandomVariate[NormalDistribution[0.,15.]]},{i,1,20}] ;
{xdat2,ydat2}=Transpose[noisydata] ;
datplot2P=Graphics[{Hue[0.65,1,1.],PointSize[0.03],Point /@ noisydata}];
datplot2=Show[ListLinePlot[noisydata, PlotStyle-> Hue[0.6,1.,1.], PlotRange-> All], datplot2P];
Show[datplot2,datplot1]

Comparing the two data (blue vs. red), we see that one is noisier than the other. So if we try to find a relationship between x and y, which one wouild be more meaningful? Intuitively, we think the red data indicates a better relationship between x and y. So, in fact,  there is a measure of the "correlation" between x and y defined as:
             r=(((OverBar[x y])  - (OverBar[x])  (Overscript[y, _ ]))/((Subscript[σ, x])  (Subscript[σ, y]) ))
where (Subscript[σ, x]) and (Subscript[σ, y])  are called the standard of deviation of the data defined as:
(Subscript[σ, x]) =(Sqrt[(∑_(i=1)^n)((((Subscript[x, i]) -(OverBar[x]) ))^2)/n ])     and (Subscript[σ, y])  =(Sqrt[(∑_(i=1)^n)((((Subscript[y, i]) -(OverBar[y]) ))^2)/n ])    
(note: this is different from unbiased SD, which is divided by (n-1)).
The closer r is to 1, the stronger the correlation between the variable.

{μX,μY}=Mean /@ ({xD,yD}/.  { xD-> xdat1 , yD->  ydat1});
{σX,σY}={(Sqrt[Mean[(xD-μX)^2]])/. xD-> xdat1 ,(Sqrt[Mean[(yD-μY)^2]])/. yD->  ydat1};
r1=(((Mean[xD*yD] - Mean[xD]*Mean[yD]))/(σX*σY))/.  { xD-> xdat1 , yD->  ydat1};
Grid[{{Style["correlation coeff r for data1=  ",24,Red,FontFamily->"Calibri"],r1}},BaseStyle->Red]

We see that the correlation is very close to 1, indicating a strong correlation of x and y. This means that the linear fit we had is a good one.
Let's do the same thing for the other data

{μX,μY}=Mean /@ ({xD,yD}/.  { xD-> xdat2 , yD->  ydat2});
{σX,σY}={(Sqrt[Mean[(xD-μX)^2]])/. xD-> xdat2,(Sqrt[Mean[(yD-μY)^2]])/. yD->  ydat2};
r2=(((Mean[xD*yD] - Mean[xD]*Mean[yD]))/(σX*σY))/.  { xD-> xdat2, yD->  ydat2};
Grid[{{Style["correlation coeff r for data2=  ",24,Blue,FontFamily->"Calibri"],r2}},BaseStyle->Blue]

Show[datplot2,datplot1
,Epilog-> {Inset["r1="<>ToString[r1],Scaled[{0.2,0.9}],BaseStyle->{20, Red}]
,Inset["r2="<>ToString[r2],Scaled[{0.2,0.8}],BaseStyle->{20, Blue}]}]

We see that the correlation of blue data is not as good as the previous red data, indicating a weaker correlation between x and y. The effect of noise (or interference) on the correlation coefficient will be studied in HW.

2. General least-square fit

Assignments and demo in class: use the app to study:
- Regression statistics
- Residue
- Coefficient confidence
- ANOVA (analysis of variance)

3. Linearization of nonlinear data

3.1 Basic idea

What if x and y are not really linear dependent of each other. For example,we have the following plot of the NASDAQ

Nasdaq={755.2`,793.73`,817.21`,843.98`,864.58`,933.45`,1001.21`,1020.11`,1043.54`,1036.06`,1059.2`,1052.13`,1059.79`,1100.05`,1101.4`,1190.52`,1243.43`,1185.02`,1080.59`,1141.5`,1226.92`,1221.51`,1292.61`,1291.03`,1379.85`,1309.`,1221.7`,1260.76`,1400.32`,1442.07`,1593.81`,1587.32`,1685.69`,1593.61`,1600.55`,1570.35`,1619.36`,1770.51`,1835.68`,1868.41`,1778.87`,1894.74`,1872.39`,1499.25`,1693.84`,1771.39`,1949.54`,2192.69`,2505.89`,2288.03`,2461.4`,2542.86`,2470.52`,2686.12`,2638.49`,2739.35`,2746.16`,2966.43`,3336.16`,4069.31`,3940.35`,4696.69`,4572.83`} ;
ListLinePlot[Nasdaq]

As a function of time, it grows faster than linear, so let's try log plot

logNAQ=Log[Nasdaq];
Nasdaqplot=ListLinePlot[logNAQ,Frame-> True, GridLines->Automatic]

That is a bit more reasonable to fit a straightline to. So let's do a linear regression with log(y)

3.2 Example of log fit

logNAQ=Log[Nasdaq];
xdat=Range[Length[Nasdaq]] ;
a=((Mean[xdat*logNAQ] -Mean[xdat]*Mean[logNAQ])/(Mean[xdat^2]-Mean[xdat]^2));
b=Mean[logNAQ] - a*Mean[xdat];
Grid[Join[{{Style["Linear regression coef",Bold],"estimate"}}
,MapThread[Prepend,{{{a},{b}}
,{"a (slope)","b (intercept)"}}]]
,Frame-> All,BaseStyle->{20,FontFamily-> "Calibri"}
,Background->{{Hue[0.16,1,1],},{Hue[0.4,1,1],}}]

NAQfit=Plot[a*x+b,{x,1,Length[xdat]}, PlotStyle->{RGBColor[1.,0.,0]}
,Frame-> True, GridLines->Automatic] ;
Show[Nasdaqplot,NAQfit]

What does this tell us:
we have:
Log(Nasdaq) =  0.0234*time (monthly from Jan 1995)) + 6.623
Therefore:
Nasdaq = (e^(0.0234*t+6.623)) = 752 *  (e^(0.0234*t)) 
It means that the Nasdaq was growing at a rate of 2.34% monthly average from Jan 1995 to March 2000, (before it tumbled down). Is 32% annual growth rate sustainable for 5 years?

3.3 Others

In general, x and y can have other nonlinear relationships, but if we know it, we can transform both x and y into variables that have linear relationship with each other, then we can perform linear regression.
 See example in HW.

4. Multivariate linear regression

4.1 Basic idea

What if y does not depend on just one parameter but many? Such as:
y = (Subscript[a, 0]) +  (Subscript[a, 1]) (Subscript[x, 1]) +  (Subscript[a, 2]) (Subscript[x, 2]).
The same principle of minimization of the error square can be applied to obtain a least square. Here, the solution for fitting parameters a's can be put in a matrix form:
 (({{1, (OverBar[(Subscript[x, 1])]), (OverBar[(Subscript[x, 2])])}, {(OverBar[(Subscript[x, 1])]), (OverBar[((Subscript[x, 1])^2)]), (OverBar[(Subscript[x, 1])(Subscript[x, 2])])}, {(OverBar[(Subscript[x, 2])]), (OverBar[(Subscript[x, 1])(Subscript[x, 2])]), (OverBar[((Subscript[x, 2])^2)])}}))(({{(Subscript[(OverHat[a]), 0])}, {(Subscript[(OverHat[a]), 1])}, {(Subscript[(OverHat[a]), 2])}})) =  (({{(OverBar[y])}, {(OverBar[(Subscript[x, 1])y])}, {(OverBar[(Subscript[x, 2])y])}})) 
 (({{(Subscript[(OverHat[a]), 0])}, {(Subscript[(OverHat[a]), 1])}, {(Subscript[(OverHat[a]), 2])}}))=(((({{1, (OverBar[(Subscript[x, 1])]), (OverBar[(Subscript[x, 2])])}, {(OverBar[(Subscript[x, 1])]), (OverBar[((Subscript[x, 1])^2)]), (OverBar[(Subscript[x, 1])(Subscript[x, 2])])}, {(OverBar[(Subscript[x, 2])]), (OverBar[(Subscript[x, 1])(Subscript[x, 2])]), (OverBar[((Subscript[x, 2])^2)])}})))^(-1)). (({{(OverBar[y])}, {(OverBar[(Subscript[x, 1])y])}, {(OverBar[(Subscript[x, 2])y])}}))

4.2 Example 1- numerical

4.2.1 illustration 1

Manipulate[
multiydat=Table[ 2.2+5*i-1.1*j+RandomVariate[LaplaceDistribution[0.,σnoise]],{i,1,10},{j,2,20,2}] ;
ListPlot3D[multiydat, BoxRatios->{1,1,1}]
,{{σnoise,2.},0.5,5}]

4.2.2 illustration 2: varying model parameters

Manipulate[
If[TrueQ[RUNmultivar1],
Column[{
Manipulate[
If[!TrueQ[{a0old,a1old,a2old,σnoiseold}=={a0,a1,a2,σnoise}],
x1arr=Sort[RandomReal[{0.,1.},10]];Dx1=Max[x1arr]-Min[x1arr];
x2arr=Sort[RandomReal[{0.,2.},15]];Dx2=Max[x2arr]-Min[x2arr];
yarr=Table[ a0+a1   x1arr[[i]]+a2*x2arr[[j]]
 +RandomVariate[LaplaceDistribution[0.,σnoise]],{i,1,Length @ x1arr }
,{j,1,Length @ x2arr }];
ydatplot=Flatten[Table[{x1arr[[i]],x2arr[[j]],yarr[[i,j]]},{i,1,Length @ x1arr}
,{j,1,Length @ x2arr}],1];
{x1R,x2R,yR}=Transpose[ydatplot];
Dy=Max[yR]-Min[yR];
matR=(({{1, Mean[x1R], Mean[x2R]}, {Mean[x1R], Mean[x1R^2], Mean[x1R*x2R]}, {Mean[x2R], Mean[x1R*x2R], Mean[x2R^2]}}));
matRinv=Inverse[matR]; Uvec=(({{Mean[yR]}, {Mean[x1R*yR]}, {Mean[x2R*yR]}}));
(({{a0E}, {a1E}, {a2E}})) = matRinv . Uvec;
{a0old,a1old,a2old,σnoiseold}={a0,a1,a2,σnoise}
,];

Row[{Column[{ ListPlot3D[ydatplot,  BoxRatios->{1,1,0.65},ImageSize->350]
,Show[Plot3D[a0E+a1E u1+a2E u2,{u1,Min[x1arr],Max[x1arr]},{u2,Min[x2arr],Max[x2arr]}
,PlotRange->{{Min[x1arr]-0.1Dx1,Max[x1arr]+0.1Dx1},{Min[x2arr]-0.1Dx2,Max[x2arr]+0.1Dx2},{Min[yR]-0.1Dy,Max[yR]+0.1Dy}}
,PlotStyle-> Hue[0.5,0.75,1,0.8]]
,Graphics3D[{Hue[0.055,1.,1.,0.6],Map[Ellipsoid[#,0.022{ Dx1, Dx2,Dy/0.65}]&,ydatplot]}]
,BoxRatios->{1,1,0.65}
,ImageSize->350]
}]


,Spacer[20]
,Column[{
Grid[{{,"matrix","inverse"}
,{, TraditionalForm @ NumberForm[matR,{5,2}], TraditionalForm @  EngineeringForm[matRinv,3]}
,{,TraditionalForm @ Uvec,}}
,BaseStyle->{16,FontFamily-> "Calibri"}, Frame-> All]

,Grid[{{"set coeff","estimated val.",},
Join[TraditionalForm/@ {(({{a0}, {a1}, {a2}})),(({{a0E}, {a1E}, {a2E}}))},{}]}
,BaseStyle->{20,FontFamily-> "Calibri"}, Frame-> All]}]

}]


,Column[{Button["Model parameters",BaseStyle-> {14, Bold,Black,FontFamily-> "Arial"}
,Background-> Hue[0.16,0.7,1]]
,Grid[{{
Control[{{a0,a0init,"offset"},0.,10.,ImageSize-> Medium}]
,Control[{{σnoise,0.5,"noise"},0.02,5,ImageSize-> Medium}]}
,{Control[{{a1,a1init,"coef 1"},-5.,5.,ImageSize-> Medium}]
,Control[{{a2,a2init,"coef 2"},-5.,5.,ImageSize-> Medium}]}}]
},Frame-> True
,FrameStyle->{Thick,Hue[0.,1.,1]}
,Background-> Hue[0.1,0.1,1]]
]
}]
,Button[,{RUNmultivar1=True }]]

,Column[{
Row[{Button["ECE3340-APP-Linear regression",
BaseStyle->{18,Black, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.4,1,1]],Spacer[20]
,TextCell["Illustration of 2-variate regression",
18,Black, Italic,FontFamily-> "Calibri"
,Background->Hue[0.75,0.2,1]]}]}]
,Row[{Button["RUN STATUS/CONTROL",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUNmultivar1,False,"->"},{False,True},ControlType-> Checkbox}] 
,Spacer[35]
,Button[Dynamic[If[TrueQ[RUNmultivar1],,]]
,{ If[TrueQ[RUNmultivar1],RUNmultivar1=False,RUNmultivar1=True]}]}]
,Initialization:>{
a0init=2.; a1init=3.7; a2init=-1.2;
}
]

4.3 Example 2

See example in HW.

5. General regression of basis functions

5.1 Introduction

In a more general way, we can consider a broader view of the least square method for any function of the type:
     y=(({{(Subscript[a, 0]), (Subscript[a, 1]), (Subscript[a, 2]), ..., (Subscript[a, n])}})). (({{1}, {(Subscript[z, 1])(x)}, {(Subscript[z, 2])(x)}, {.}, {(Subscript[z, n])(x)}})) 
where (Subscript[z, 1])(x), ... (Subscript[z, n])(x) are any given functions of x, which do not have to be linear. The problem is to find the coefficients  (Subscript[a, 0]) and (({{(Subscript[a, 1]), (Subscript[a, 2]), ..., (Subscript[a, n])}})) such that we obtain a least square fit to the data.

5.2 Example 1

See example in HW.

Consider the case when we have a summing recording of three signals. Suppose we know characteristics function of each signal, (Subscript[s, i])[t], is it possible to use regression to determine the amplitude of each signal?
                   S[t] =(Subscript[a, 0])+(Subscript[a, 1]) (Subscript[s, 1])[t]+(Subscript[a, 2]) (Subscript[s, 2])[t]+(Subscript[a, 3]) (Subscript[s, 3])[t]

5.3 General formalism

The general format for least square fitting is to form the equation:
                          Y=Z ·A+e   
where Y  is the data array, Z is the matrix of the functions evaluated at each x value, A is the coefficient vector, and  e is the residual error to be minimized. The sum square of the residual error is:                 
                           (e^T). e =(((Y-Z ·A))^T). (Y-Z ·A) =(Y^T). Y -((Y^T) Z ·A +(A^T)(Z^T)Y)+(A^T)(Z^T) Z ·A

(∂/(∂(A^T)))(e^T). e =(Z^T) Z ·A-(Z^T)Y     
Minimum error requires:              (∂/(∂(A^T)))(e^T). e =0=(Z^T) Z ·A-(Z^T)Y

Or:                                           (Z^T) Z ·A=(Z^T)Y   
                                                (OverHat[A])=((((Z^T) Z ))^(-1)) (Z^T)Y

Numerically, a test value of the fit is:  (((A-(OverHat[A])) . (A-(OverHat[A])))/(A . A))

In Mathematica, we can use the Fit function:
Fit[data,{basis function list},variable]

5.4 Example 2

A function y depends on t as follow:
            y=(Subscript[a, 0]) + (Subscript[a, 1])((t^2)/(1+t)) Cos[3 π t] + (Subscript[a, 2])(e^(-t)) (t^2) (t-2)+(Subscript[a, 3]) (e^(-t/2)) Sin[8 π t]

Determine the coefficients.

Manipulate[
If[TrueQ[RUNgenfuncFit],
Column[{
Manipulate[
Avec={a0,a1,a2,a3};
If[!TrueQ[{Avecold,σnoiseold,nSold}=={Avec,σnoise,nS}] || TrueQ[NEWDAT],
tarr=Sort[RandomReal[{0.,4.},nS]];tR=Max[tarr]-Min[tarr];
ZmatT={ConstantArray[1,nS],fB1[#] &/@ tarr,fB2[#] &/@ tarr,fB3[#] &/@ tarr};
Zmat=Transpose[ZmatT];
yarrF=Table[ Avec . Zmat[[i]]
 +RandomVariate[NormalDistribution[0.,σnoise]],{i,nS }];
datplot=Transpose[{tarr,yarrF}];
Aest=Inverse[ZmatT. Zmat]. ZmatT. yarrF;
ρE=(((Avec-Aest) . (Avec-Aest))/(Avec . Avec)) ;
outfit=Fit[datplot,{1,fB1[t],fB2[t],fB3[t]},t];
{Avecold,σnoiseold,nSold}={Avec,σnoise,nS};NEWDAT=False
,];
Column[{
Row[{
Show[Plot[a0+a1 fB1[t]+a2 fB2[t]+a3 fB3[t],{t,0,4}
,PlotRange-> {{0,4},All}
,Frame-> True, GridLines-> Automatic
,PlotStyle->{Hue[0.8,1.,0.75]}
,Filling-> {1-> {0.,Hue[0.77,0.7,1,0.4]}}]
,ListLinePlot[datplot,PlotStyle-> {Thickness[0.001],Black}]
,Graphics[{Red,PointSize[0.02],Point /@ datplot }]
,fBplot,PlotRange-> {{0,4},All}
,Epilog-> Map[Inset[TraditionalForm @ #[[1]],#[[2]],BaseStyle->{20,Bold}]&
,{{"(Subscript[z, 1])(t)",{3.7,-6.3}},{"(Subscript[z, 2])(t)",{3.7,-4.5}},{"(Subscript[z, 3])(t)",{3.7,-2.5}}}]
,AspectRatio->1.6,ImageSize-> 300]
,Spacer[20]
,Column[{

,Grid[{{"set coeff","estimated val.",},
Join[TraditionalForm/@ {(({{a0}, {a1}, {a2}, {a3}})),Table[{Aest[[m]]},{m,1,Length @ Aest}]},{}]
,{, ρE,}}
,BaseStyle->{18,FontFamily-> "Calibri"}, Frame-> All]
}]

}]
,Grid[{{TextCell["Mathematica Fit function output:",18,Italic]}
,{Style[TraditionalForm @ outfit,18,FontFamily->"Times New Roman",Bold,Background->Hue[0.12,0.5,1]]}}]
}]

,Column[{Button["Model parameters",BaseStyle-> {14, Bold,Black,FontFamily-> "Arial"}
,Background-> Hue[0.16,0.7,1]]
,Grid[{{
Control[{{a0,a0init,"offset"},0.,10.,ImageSize-> Medium}]
,Control[{{σnoise,0.5,"noise"},0.02,5,ImageSize-> Medium}]
,Row[{Control[{{nS,25,"#sampl. pnts"},{10,25,50,100},ControlType-> SetterBar}]
,Spacer[5],Button[,NEWDAT=True]}]}
,{Control[{{a1,a1init,"coef 1"},-5.,5.,ImageSize-> Medium}]
,Control[{{a2,a2init,"coef 2"},-5.,5.,ImageSize-> Medium}]
,Control[{{a3,a3init,"coef 3"},-5.,5.,ImageSize-> Medium}]}}]
},Frame-> True
,FrameStyle->{Thick,Hue[0.,1.,1]}
,Background-> Hue[0.1,0.1,1]]
]
}]
,Button[,{RUNgenfuncFit=True }]]

,Column[{
Row[{Button["ECE3340-APP-Linear fit with basis functions",
BaseStyle->{18,Black, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.4,1,1]],Spacer[20]
,TextCell["Illustration of Fit function",
18,Black, Italic,FontFamily-> "Calibri"
,Background->Hue[0.75,0.2,1]]}]}]
,Row[{Button["RUN STATUS/CONTROL",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUNgenfuncFit,False,"->"},{False,True},ControlType-> Checkbox}] 
,Spacer[35]
,Button[Dynamic[If[TrueQ[RUNgenfuncFit],,]]
,{ If[TrueQ[RUNgenfuncFit],RUNgenfuncFit=False,RUNgenfuncFit=True]}]}]
,Initialization:>{
a0init=1.5; a1init=1.; a2init=4;a3init=-3.;
fB1[t_]:=((t^2)/(1+t)) Cos[3 π t]; fB2[t_]:=(E^(-t)) (t^2) (t-2); fB3[t_]:=(E^(-t/2)) Sin[8 π t] ;
fBplot=Plot[{fB1[t]-8,fB2[t]-5,fB3[t]-3},{t,0,4}
,PlotStyle->{Hue[0.05,1.,1],Hue[0.35,1.,0.6],Hue[0.6,1.,1]}
,Filling-> {1-> {-8,Hue[0.05,0.8,1,0.4]}
,2-> {-5,Hue[0.35,0.8,1,0.4]},3-> {-3,Hue[0.6,0.8,1,0.4]}}
,Epilog-> Map[Inset[TraditionalForm @ #[[1]],#[[2]]]&
,{{(Subscript[z, 1])(x),{3.5,-6}},{(Subscript[z, 2])(x),{3.5,-5}},{(Subscript[z, 3])(x),{3.5,-3}}}]]
}
]

6. General (non-linear) model parameter fit

6.1 Discussion

Consider a general function:     f[x ; {(Subscript[a, 0]),(Subscript[a, 1]),..., (Subscript[a, p])}]  
where          x={(Subscript[x, 1]), (Subscript[x, 2]), ... (Subscript[x, n])}
and              {(Subscript[a, 0]),(Subscript[a, 1]),..., (Subscript[a, p])} are unknown parameters
Suppose we have data:
           (Subscript[y, k])->f[(Subscript[x, k]) ; {(Subscript[a, 0]),(Subscript[a, 1]),..., (Subscript[a, p])}]   for k from 1 to m
Then, parameters   {(Subscript[a, 0]),(Subscript[a, 1]),..., (Subscript[a, p])} can be estimated by mininize a metric of          M[{(Subscript[y, k])} , f[(Subscript[x, k]) ; {(Subscript[a, 0]),(Subscript[a, 1]),..., (Subscript[a, p])}] ] 
Usually, least-square is used, and estimates for parameters are obtain as a LQ-minimization problem:
       M[{(Subscript[a, 0]),(Subscript[a, 1]),..., (Subscript[a, p])}]=(∑_(k=1)^m)((((Subscript[y, k])-f[(Subscript[x, k]) ; {(Subscript[a, 0]),(Subscript[a, 1]),..., (Subscript[a, p])}]))^2)

Linear regression is just a special case of this, in which the model function is linear and can be solved with linear algebra methods.

In principle, the minimization can be done with a variety algorithms, most common is the method of steepest descend.

6.2 Example: Sigmoid or logistic function

A very common function related to exponential is the sigmoid or logistic function that appears often in both physical science and economics is the sigmoid or logistic function (also known as Fermi-Dirac in physics).

(Subscript[f, FD])[t,s,(Subscript[T, 0]),β]=(s/(1+(e^(-(t-(Subscript[T, 0]))/β)))) 
                                  (Subscript[f, cFD])[t,s,(Subscript[T, 0]),β]=(s/(1+(e^((t-(Subscript[T, 0]))/β))))

Manipulate[
Plot[{(1/(1+(E^(-x/β)))),(1/(1+(E^(x/β))))},{x,-8,8}]
,{{β,1},0.1,2}]

This function is the essence behind a popular regression method for categorical variables, known as logistic regression.

6.2.1 Technology adopter model

Below is the data of smartphone quarterly shipment as a function of time. What is the annual growth rate of smart phone shipment over a 4-yr period from 2009 to 2012?  We can answer the question by plotting the data on a semilog scale from Q1:09 to Q4:12 (4-year period). From the log scale, as above, we can easily determine the annual adoption growth rate.

Another example:

However, can these be exponential forever? Of course not. Ultimately, it will saturate at 100% of the population. This is known as early adoption exponential, linear growth, and saturation slow down.

Early adoption growth tends to be exponential, followed by a brief constant-like growth, and the end-game trend is saturation (when everyone already has enough!). The behavior is very common throughout human history and the curve to describe it is call sigmoid or logistic curve (S-curve, etc.) , which is mathematically very similar to or of the same form as Ferm-Dirac distribution (reverse for growth)

Source code

Manipulate[
If[TrueQ[LOGISTICDEMO==]
,
Manipulate[Tm=τ Log[(E^(T/τ)) T] ;
{Show[Plot[ 100.t ((E^((t-Tm)/τ))/(1+ t (E^((t-Tm)/τ)))),{t,0,20.}, PlotRange->{0.,100.}, Frame-> True
,GridLines-> Automatic, PlotStyle-> {Thick, Red}, Filling-> {1-> {Axis, Hue[0.3,1,1,0.5]}}
,FrameLabel-> {"Time","% of adoption"}
,LabelStyle->{16,FontFamily-> "Arial"}
,Epilog->{ Inset[Style["Early\nadopters",20, Red
,FontFamily-> "Calibri",LineSpacing->{0,20}]
,ImageScaled[{0.25,0.35}]]
,Inset[Style["Main\nstreamers",20, Hue[0.35,1,0.5]
,FontFamily-> "Calibri",LineSpacing->{0,20}],ImageScaled[{0.50,0.55}]]
, Inset[Style["Late\nadopters",20,Blue
,FontFamily-> "Calibri",LineSpacing->{0,20}],ImageScaled[{0.850,0.85}]]}]
,Graphics[{Hue[0.8,1,1],Dashed, Thickness[0.005], Line[{{T,0.},{T,100.}}]}]
,ImageSize-> {400,270}]
,Show[LogPlot[100. t ((E^((t-Tm)/τ))/(1+ t (E^((t-Tm)/τ)))),{t,0.15 T,20.}, PlotRange->{{0.,20.},{(10^(-1)),100.}}, Frame-> True
,GridLines-> Automatic, PlotStyle-> {Thick, Red}, Filling-> {1-> {Axis, Hue[0.3,1,1,0.5]}}
,FrameLabel-> {"Time","% of adoption"}
,LabelStyle->{16,FontFamily-> "Arial"}
,Epilog->{Inset[Style["Early\nadopters",20, Red
,FontFamily-> "Calibri",LineSpacing->{0,20}]
,ImageScaled[{0.25,0.3}]] 
, Inset[Style["Late\nadopters",20,Blue
,FontFamily-> "Calibri",LineSpacing->{0,20}]
,ImageScaled[{0.80,0.85}]]}]
,Graphics[{Hue[0.8,1,1],Dashed, Thickness[0.005], Line[{{T,Log[0.01]},{T,Log[100.]}}]}]
,ImageSize-> {400,270}]
}
,{{τ,0.75, "Adoption growth time"},0.1,2.5} 
,{{T,10., "Mean time to adoption"},2.,20.}  ]
,]
,Grid[{{Control[{{LOGISTICDEMO,,},{,},ControlType->Toggler}] 
,Item[Style["Sigmoid/Logistic/S-curve demo of technology adoption",16,FontFamily->"Century Gothic"]
,Background->Hue[0.75,0.35,1]]
}}]]

Demo

Fermi-Dirac distribution for comparison

6.2.2 Technology phase-out model

The complementary from of upward S-curve is the downward S-curve which is:
                             1-(1/(1+(e^(-(t-(Subscript[T, 0]))/β))))=(1/(1+(e^((t-(Subscript[T, 0]))/β))))

6.2.3 See homework

Pick 2 technologies in APP Sigmoid model fit for technology adoption of American households (below),fit with the S-curve to obtain mean time to adoption and the adoption annual growth rate.