>	read("/home/scott/www/mat331.spr08/problems/lsq_data.txt");

 

defined line_pts(), bad_line_pts(), quadratic_pts(), cubic_pts(), and circle_pts()

 

>	stuff:=bad_line_pts():

 

>	plot(stuff,style=point);

 

>	badline:=CurveFitting[LeastSquares](stuff,x);

 

(1)

 

>	plot([stuff,badline],x=-10..10,style=[point,line]);

 

>	nops(stuff);

 

(2)

 

>	epsilon:=(m,b,pt) -> ( m*pt[1]+b - pt[2]);

 

(3)

 

>	H:=(m,b,data) -> sum( epsilon(m,b,data[i])^2,i=1..nops(data));

 

(4)

 

>	H(1,2,stuff);

 

(5)

 

>	plot3d(H(m,b,stuff),m=-3..0,b=-10..0, axes=boxed);

 

Rachel suggests 1-exp(-x^2).  Let's try it. 

>	R:=(m,b,data) -> sum( 1-exp(-epsilon(m,b,data[i])^2),i=1..nops(data));

 

(6)

 

>	plot3d(R(m,b,stuff),m=-3..0,b=-10..0, axes=boxed);

 

>	plot3d(R(m,b,stuff),m=-3..-1,b=-10..-7.5, axes=boxed);

 

>	solve({diff(R(m,b,stuff),m)=0, diff(R(m,b,stuff),b)=0}, {m,b});

 

Warning,  computation interrupted

 

>	fsolve({diff(R(m,b,stuff),m)=0, diff(R(m,b,stuff),b)=0}, {m,b}, {m=-2.5..-1.5, b=-10..-9});

 

(7)

 

>	Rline:=subs(%,m*x+b);

 

(8)

 

>	plot([stuff,badline,Rline],x=-10..10,style=[point,line,line]);

 

>	plot([stuff,badline,Rline],x=-10..10,y=-40..20,style=[point,line,line]);

 

>	S:=(m,b,data) -> sum( ln(1+epsilon(m,b,data[i])^2),i=1..nops(data));

 

(9)

 

>	plot3d(S(m,b,stuff),m=-3..0,b=-20..0, axes=boxed);

 

>	plot3d(S(m,b,stuff),m=-2.2..-2,b=-10..-9, axes=boxed);

 

>	fsolve({diff(S(m,b,stuff),m)=0, diff(S(m,b,stuff),b)=0}, {m,b}, {m=-2.5..-1.5, b=-10..-9});

 

(10)

 

>	Sline:=subs(%,m*x+b);

 

(11)

 

>	plot([stuff,badline,Rline,Sline],x=-10..10,y=-40..20,style=[point,line,line,line]);

 

>

?Minimize

 

>	Optimization[Minimize](S(m,b,stuff));

 

(12)

 

>	Optimization[Minimize](R(m,b,stuff));

 

(13)

 

>	T:=(m,b,data) -> sum( abs(epsilon(m,b,data[i])),i=1..nops(data));

 

(14)

 

>	T(1,2,stuff);

 

(15)

 

>	plot3d(T(m,b,stuff),m=-3..0,b=-20..0, axes=boxed);

 

>	Optimization[Minimize](T(m,b,stuff));

 

Error, (in Optimization:-NLPSolve) no improved point could be found

 

>	plot3d(T(m,b,stuff),m=-2.25..-2,b=-9.5..-9.25, axes=boxed);

 

>

 

>