> points := []:
> points := [op(points), [40, 368.0]]:
> points := [op(points), [100, 2447.0]]:
> points := [op(points), [200, 10051.0]]:
> points := [op(points), [300, 23846.0]]:
> points := [op(points), [500, 64112.0]]:
> points := [op(points), [2000, 960156.0]]:
> points;
> 

  [[40, 368.0], [100, 2447.0], [200, 10051.0], [300, 23846.0],

        [500, 64112.0], [2000, 960156.0]]

> m := nops(points);

                                m := 6

> x_coords := []: y_coords:=[]: yw_coords:=[]:
> 
> for i from 1 to m do
>     x_coords := [op(x_coords), points[i,1]];
>     y_coords := [op(y_coords), points[i,2]];yw_coords :=
> [op(yw_coords), Weight(points[i,2],1/points[i,1]^2)];
> od:
> x_coords;

                    [40, 100, 200, 300, 500, 2000]

> y_coords;

         [368.0, 2447.0, 10051.0, 23846.0, 64112.0, 960156.0]

> model := stats[fit,leastsquare[[x,y],y=a*x^2 + b*x + c,
> {a,b,c}]]([x_coords, y_coords]);

                                2
      model := y = .2339439835 x  + 12.72771129 x - 1062.235228

> with(linalg):
> A := matrix(m,3):
> for i from 1 to m do
>    for j from 0 to 2 do A[i,j+1] := x_coords[i]^j;
> od: od:
> print(A);

                        [1      40       1600]
                        [                    ]
                        [1     100      10000]
                        [                    ]
                        [1     200      40000]
                        [                    ]
                        [1     300      90000]
                        [                    ]
                        [1     500     250000]
                        [                    ]
                        [1    2000    4000000]

> linsolve( evalm(transpose(A) &* A), evalm(transpose(A) &*
> vector(y_coords)) );

               [-1062.235432, 12.72771283, .2339439827]

> floor(subs(x=10000,op(2,model)));

                               23520613

> yw_coords;

  [Weight(368.0, 1/1600), Weight(2447.0, 1/10000),

        Weight(10051.0, 1/40000), Weight(23846.0, 1/90000),

        Weight(64112.0, 1/250000), Weight(960156.0, 1/4000000)]

> modelw := stats[fit,leastsquare[[x,y],y=a*x^2 + b*x + c,
> {a,b,c}]]([x_coords, yw_coords]);

                                 2
      modelw := y = .2367384129 x  + 7.478153418 x - 360.3026286

> floor(subs(x=10000,op(2,modelw)));

                               23748262

>