I am trying a sample optimization problem. I create 150 points on an
ellipse in the x-y plane. The major axis of the ellipse is x, and the
minor is y (the two foci of the ellipse lie on the x-axis at equal
distances from the origin). Then, I rotate the axes by -30°, and then
shift the points by (0.25,0.25) and introduce some noise. After that,
I try to recover the ellipse from the noisy data.
The model I use to try to recover the ellipse from the noisy data is
based on the fact that if the ellipse is shifted back to the origin,
the eigenvectors of the resulting quadratic will each lie along a
major axis. For the ellipse which is shifted to the origin (but still
tilted) the equation is w'Aw=1, where A is the matrix describing the
quadratic form, w is a vector from the origin to any point on the
ellipse, and the prime denotes "transpose" (as in Matlab/Scilab).
Since the 2×2 matrix A is symmetric, it is diagonalizable, has two
distinct (and orthogonal) eigenvectors e1 and e2, with corresponding
eigenvalues λ1 and λ2. Therefore, A=λ1*e1*e1' + λ2*e2*e2' (observe
that the ROW eigenvector is the SECOND vector in each product; for
details see Gilbert Strang's book "Linear Algebra and its
Applications" and/or watch his online lectures at ocw.mit.edu). The
model is based on this.
Here is the code:
import numpy,scipy,os,sys,pylab,matplotlib
from scipy import pi as PI
# Create 150 points on the ellipse x**2+4*y**2=1
phi=scipy.linspace(-PI,+PI,150)
theta=var('theta')
theta=PI/6.
x=numpy.cos(phi);y=0.5*numpy.sin(phi)
# Rotate the axes by -30° 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)]
# Try to recover the tilted and shifted ellipse from the noisy data
x1,y1,x2,y2,lambda1,lambda2,X,Y=var
('x_1,y_1,x_2,y_2,lambda_1,lambda_2,X,Y')
B=transpose( numpy.array([xx,yy,numpy.ones
(150)],dtype='float64') )
X=1.0
Y=1.0
lambda1=1.0
lambda2=1.0
x1=1.0
y1=1.0
x2=1.0
y2=1.0
xx=xx.reshape(150,1);yy=yy.reshape(150,1)
def elp(xx,yy,X,Y,x1,y1,x2,y2,lambda1,lambda2):
xo=xx-numpy.ones((150,1))*X #x-offsets
yo=yy-numpy.ones((150,1))*Y # y-offsets
return lambda1*((x1*xo)**2+(y1*yo)**2) + lambda2*( (x2*xo)**2+
(y2*yo)**2 )
find_fit(B,elp,parameters=[X,Y,x1,y1,x2,y2,lambda1,lambda2],variables=
[xx,yy],solution_dict=1)
When I run the above code, I get:
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/home/nonsense/.sage/sage_notebook/worksheets/linuxgus/1/code/
1.py", line 38, in <module>
find_fit(B,elp,parameters=
[X,Y,x1,y1,x2,y2,lambda1,lambda2],variables=
[xx,yy],solution_dict=_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/sage/
numerical/optimize.py", line 553, in find_fit
estimated_params, d = leastsq(error_function, initial_guess, args
= (x_data, y_data))
File "/usr/local/sage-3.3/local/lib/python2.5/site-packages/scipy/
optimize/minpack.py", line 264, in leastsq
m = check_func(func,x0,args,n)[0]
File "/usr/local/sage-3.3/local/lib/python2.5/site-packages/scipy/
optimize/minpack.py", line 11, in check_func
res = atleast_1d(thefunc(*((x0[:numinputs],)+args)))
File "/usr/local/sage-3.3/local/lib/python2.5/site-packages/sage/
numerical/optimize.py", line 546, in error_function
result[row] = func(*fparams)
ValueError: setting an array element with a sequence.
Can somebody shed some light on this and/or suggest alternative ways
of doing this?
--~--~---------~--~----~------------~-------~--~----~
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
-~----------~----~----~----~------~----~------~--~---
