Problems (?) with least squares

[linuxgus]

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)

p0=[1.r,1.r,1.r,1.r,1.r,1.r] #initial guess
leastsq(residual,p0,args=(xx,yy),full_output=1)


(In the above code I omitted some plotting commands for clarity.)

I get the following error (expanded):
ValueError: setting an array element with a sequence.
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/home/gus/.sage/sage_notebook/worksheets/linuxgus/1/code/
86.py", line 16, in <module>
leastsq(residual,p0,args=(xx,yy),full_output=_sage_const_1 )
File "/usr/local/sage-3.3/local/lib/python2.5/site-packages/
SQLAlchemy-0.4.6-py2.5.egg/", line 1, in <module>

File "/usr/local/sage-3.3/local/lib/python2.5/site-packages/scipy/
optimize/minpack.py", line 268, in leastsq
retval = _minpack._lmdif
(func,x0,args,full_output,ftol,xtol,gtol,maxfev,epsfcn,factor,diag)
minpack.error: Result from function call is not a proper array of
floats.
ValueError: setting an array element with a sequence.
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/home/gus/.sage/sage_notebook/worksheets/linuxgus/1/code/
86.py", line 16, in <module>
leastsq(residual,p0,args=(xx,yy),full_output=_sage_const_1 )
File "/usr/local/sage-3.3/local/lib/python2.5/site-packages/
SQLAlchemy-0.4.6-py2.5.egg/", line 1, in <module>

File "/usr/local/sage-3.3/local/lib/python2.5/site-packages/scipy/
optimize/minpack.py", line 268, in leastsq
retval = _minpack._lmdif
(func,x0,args,full_output,ftol,xtol,gtol,maxfev,epsfcn,factor,diag)
minpack.error: Result from function call is not a proper array of
floats.

This problem is probably of my own making. How is a "proper array of
floats" defined? I am trying to recover the tilt angle, phi, the
lengths of the major and minor axes, a and b, respectively, and the x
and y offset, X and Y, respectively. I don't have to use
scipy.optimize for this. If cvxopt does the job, I can use cvxopt or
any other package.

I am using sage 3.4 notebook on SuSE 11.1, 64-bit, with kernel
2.6.28.2, if that may help.

Thank you in advance for your response.

Gus

[Jason Grout]

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

[linuxgus]

On Mar 30, 7:49 am, Jason Grout <jason-s...@creativetrax.com> wrote:
> linuxgus wrote:

...

> >     """
> >     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

--Hi, Jason. Thank you for responding.

I tried your suggestion at work, where I am running Sage 3.3 (I will
probably upgrade to 3.4 today). I get
TypeError: data type not understood

I used 'float' (in single quotes) and got the same error. It is
certainly different than the error I got on Sage 3.4 at home. What is
odd, if I run the same code (without the 'type=float') on sage 3.3 at
work, I still get the same "TypeError: data type not understood", so
Sage 3.3 and 3.4 behave differently in this respect (shouldn't both
versions give me the same error?).

Gus