linuxgus wrote:
> I am trying to use least squares in sage. I started with a simple
> example: Take an ellipse, whose major and minor axes are x and y,
> respectively, turn in by -30 degrees, shift it, inject some noise, and
> then try to recover it via least squares. All examples I have seen
> deal with a single variable (noisy y data vs. x). I tried leastsq
> from scipy like this (this is an excerpt from a larger program in
> notebook; not all variables are used by this excerpt):
>
>
> import numpy,scipy,os,sys,pylab,matplotlib
> from scipy import pi as PI
> from scipy.optimize import leastsq
>
> # Create the original ellipse
> phi=scipy.linspace(-PI,+PI,150)
> theta=var('theta')
> theta=PI/6.0
> x=numpy.cos(phi);y=0.5*numpy.sin(phi)
> theta
>
>
> # Rotate the axes by -30deg and shift them by (-0.25,-0.25); inject
> some noise
> xx=x*cos(theta)+y*sin(theta)+[normalvariate(0.25,.0075) for i in range
> (150)]
> yy=-x*sin(theta)+y*cos(theta)+[normalvariate(0.25,.0075) for i in range
> (150)]
> # The xx and yy arrays hold the x and y values of the noisy, tilted,
> and shifted ellipse.
>
> # Try to recover the tilted and shifted ellipse from the noisy data
> x1,y1,x2,y2,l1,l2,X,Y,u,v=var('x1,y1,x2,y2,l1,l2,X,Y,u,v')
>
> def residual(x,y,p):
> """ X & Y are offsets
> """
> a=p[0];b=p[1];phi=p[2];delta=p[3];X=p[4];Y=p[5]
> return numpy.array(x-a*numpy.cos(phi+delta)+X , y-b*numpy.sin(phi
> +delta)+Y)
Here's a guess (I haven't checked it):
make the above function return:
return numpy.array(x-a*numpy.cos(phi+delta)+X ,
y-b*numpy.sin(phi+delta)+Y, type=float)
Jason
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/sage-support
URLs: http://www.sagemath.org
-~----------~----~----~----~------~----~------~--~---
