SimpleFit[ ] Command
********************
A modified version of the Mathematica Fit[ ] routine for doing a least-square fit of a set of basis functions to a list of data values. This modified version, called SimpleFit[], is presented here with a short example of its application to a hermite-Gaussian-squared mode expansion.
Let's define the set of data points which we want to match to be on a
broadened unit-area gaussian function with gaussian spot size w0=2
between x=0 and x=3:
*
* IN[ ]:
*
* w0=2; g[x_] := Sqrt[2./(Pi w0^2)] Exp[-2 x^2/w0^2];
*
* hgData = Table[{x, N[g[x]]}, {x, 0, 3, 0.3}];
*
* OUT[ ]:
*
* {{0, 0.399}, {0.3, 0.381}, {0.6, 0.333}, {0.9, 0.266},
* {1.2, 0.194}, {1.5, 0.13}, {1.8, 0.079}, {2.1, 0.044},
* {2.4, 0.0224}, {2.7, 0.0104}, {3., 0.00443}}
*
We want to match this data to a basis set of Hermite-gaussian squared
functions h[n,x] in the form
c[0] h[0,x] + c[1] h[1,x] + c[2] h[2,x]
where the basis functions h[n,x] are normalized Hermite-gaussian (HG)
squared intensity profiles with spot size w = Sqrt[2]. (These
correspond to the intensity profiles of lowest and higher order modes
in a stable-cavity laser device.) Note that the HG functions
themselves are orthogonal, but the HG-squared functions h[n,x] are
not. These functions as defined as follows:
*
* IN[ ]:
*
* dum = Sqrt[1./N[Pi]];
*
* h[n_,x_] = (dum/(2^n n!)) HermiteH[n,x]^2 Exp[-x^2]
*
* OUT[ ]:
* 2
* 0.56419 HermiteH[n, x]
* -----------------------
* 2
* n x
* 2 E n!
*
One way to do this expansion would be to use the Fit[] command from
Mathematica as follows. We'll use just three basis functions in this
example, and call the result hgStandardFit[x]:
*
* IN[ ]:
*
* hgBasis = Table[h[n,x], {n, 0, 2}];
*
* hgStandardFit[x_] = Fit[hgData, hgBasis, x]
*
* OUT[ ]:
*
* 2 2 2
* 0.376547 0.257976 x 0.00631176 (-2 + 4 x )
* -------- + ----------- + -----------------------
* 2 2 2
* x x x
* E E E
*
The problem with this result is that the expansion is expressed
not in the normalized functions we want, but rather in terms of
the underlying algebraic form of the Hermite functions. As a
result the numerical coefficients are not the c[n] values that
we want.
Withoff has defined as an alternative the SimpleFit[] command:
*
* IN[ ]:
*
* SimpleFit[data_, basis_, vars_] :=
* Block[{purebasis, designmatrix, designinverse, response,
* fitcoefficients},
* purebasis = Function[Evaluate[vars], Evaluate[basis]];
* designmatrix = Apply[purebasis, data, 1];
* designinverse = PseudoInverse[designmatrix];
* response = Last /@ data;
* fitcoefficients = designinverse . response;
* fitcoefficients
* ]
*
Applying this command to the same data and basis set, produces a list
hgCoefficients which contains the desired expansion coefficients we
want:
*
* IN[ ]:
*
* hgCoefficients = SimpleFit[hgData, hgBasis, x]
*
* OUT;
*
* {0.667412, 0.228625, 0.0894984}
*
The list hgCoefficients contains exactly the coefficients
c[n] that we want. In simple cases like this one, applying a
normalized basis set to a normalized (i.e., unity area) input
function as above should also yield a "total power" given by the
sum of the c[n] coefficients which is close to unity (though
this will not always happen, especially if some of the c[n] turn
out to be negative). Let's check:
*
* IN[ ]:
*
* pwr = Apply[Plus,hgCoefficients]
*
* OUT[ ]:
*
* 0.985536
*
This is not bad, and gets closer to unity if you use more
functions in the basis set.
The fitting function (call it hgSimpleFit) is then given by the
dot product of the hgCoefficients and the hgBasis, and has the
expanded form:
*
* IN[ ]:
*
* hgSimpleFit[x_] = Expand[hgCoefficients . hgBasis]
*
* OUT[ ]:
* 2 4
* 0.401794 0.156988 x 0.100988 x
* -------- + ----------- + -----------
* 2 2 2
* x x x
* E E E
*
To check that this hgSimpleFit result is the same as the
hgStandardFit result obtain using Fit[], expand the
hgStandardFit result and compare with the expanded
hgSimpleFit result:
*
* IN[ ]:
*
* Expand[hgStandardFit[x]]
*
* OUT[ ]:
*
* 2 4
* 0.401794 0.156988 x 0.100988 x
* -------- + ----------- + -----------
* 2 2 2
* x x x
* E E E
*
The expanded versions of hgStandardFit and hgSimpleFit agree,
as they should.