On Tuesday 09 December 2008 08:47:44 Valery Pipin wrote:
> Michael Carey wrote:
> > I have looked at linear regression for approximating a linear
> > graph, but how would I approximate a square or third order
> > polynomial?
> >
> > Thanks Mick C
>
> real life example. Using the Peter's Eggleton stellar evolution
> code I've got a table of the stellar interior parameters for the
> Brown Dwarf (0.2 M_sun, 0.006 L_sun). The parameters are pressure,
> temperature, luminosity, density and etc.. . Now I want to use
> their functional form in my dynamo model. The best way is to get
> the Chebyshev approximations to them. I do as follows,
> 1) read the table to maxima
> 2) find the spline interpolation
> 3) find the Chebyshev approximations
>
> The code is below. The gz-table is in attachment.
Well, here's (maybe a less practical) code with
plain polynomial approximation by weighted least squares.
--
Pungenday, Aftermath 51 YOLD 3174
Alexey Beshenov http://beshenov.ru/
-------------- next part --------------
load ("wls.mac")$
data_x : block([n:100], makelist(float(2*%pi*i/n), i,0,n))$
data_y : map('cos,data_x)$
poly_appr : least_squares(6,data_x,data_y);
plot2d(poly_appr, [x,0,2*%pi])$
-------------- next part --------------
numberlistp (l) := apply("and", map('numberp,l))$
s_pow (x,n) := if n = 0 then 1.0 else x^n$
least_squares (p, [param]) := block([weighted],
if length(param) = 2 then
weighted : false
elseif length(param) = 3 then
weighted : true
elseif not apply("and", map('listp,param)) then
error("usage: least_squares (p, [x_1, ..., x_n], [y_1, ..., y_n])
or least_squares (p, [x_1, ..., x_n], [y_1, ..., y_n], [w_1, ..., w_n])"),
if not (integerp(p) and p >= 1) then
error("p = 1, 2, ... expected"),
block([x:float(param[1]), y:float(param[2]), w, n, phi, M],
if not apply("and", map('numberlistp,param)) then
error("numeric values expected"),
if weighted then
w : float(param[3]),
n : length(y),
if length(x)#n then
error("lists of equal length expected"),
if weighted then
if length(w)#n then
error("lists of equal length expected"),
phi : zeromatrix(p+1,p+1),
M : zeromatrix(p+1,1),
for j : 0 thru p do
for k : j thru p do
if weighted then
phi[j+1][k+1] : sum(w[i]*s_pow(x[i],j)*s_pow(x[i],k), i,1,n)
else
phi[j+1][k+1] : sum(s_pow(x[i],j)*s_pow(x[i],k), i,1,n),
for j : 0 thru p do
for k : j+1 thru p do
phi[k+1][j+1] : phi[j+1][k+1],
for k : 0 thru p do
if weighted then
M[k+1][1] : sum(w[i]*y[i]*s_pow(x[i],k), i,1,n)
else
M[k+1][1] : sum(y[i]*s_pow(x[i],k), i,1,n),
b : first(linsolve_by_lu(phi,M,'floatfield)),
sum(b[i+1][1]*'x^i, i,0,p)
)
)$