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