Scipy ODR Python使用方法

Question

我正在尝试用ODR（最小二乘法的一种）把9个点拟合成一个圆锥曲线，使用的是非球面透镜的公式：

z(r) = r² /(R*(1+sqrt(1-(1+K)*(r²/R²))))

这里的R是曲率半径，K是圆锥常数，而r = sqrt(x²+y²)。K是一个固定的值（已知），而R是我想要找的值。我是从http://wiki.scipy.org/Cookbook/Least_Squares_Circle开始的，想用Python来实现这个。

我用的圆锥的隐式形式是r² - 2.R.Z + (1+K).Z²。

这是我写的代码：

# -*- coding: cp1252 -*-
from scipy import  odr
from numpy import *

# Coordinates of the 3D points
X   =   [   0,  1,  -1, 0,  0,  0.5,    -0.5,   0,  0   ]
Y   =   [   0,  0,  0,  1,  -1, 0,  0,  0.5,    -0.5    ]
Z   =   [   0,  0.113696489,    0.113696489,    0.113696489,    0.113696489,    0.027933838,    0.027933838,    0.027933838,    0.027933838]

#constantes
Rc = 8
K = -0.8

def calc_r(x, y):
    return (x**2 + y**2)

def calc_z(r, R):
    return r**2 /(R*(1+sqrt(1-(1+K)*(r**2/R**2))))

def f_3(beta, M):

    r = calc_r(M[0],M[1])
    Z = calc_z(r, beta[0])

    return r**2 - 2*beta[0]*Z + (1+K)*Z**2


beta0 = [Rc]

lsc_data  = odr.Data(row_stack([X, Y]), y=1)
lsc_model = odr.Model(f_3, implicit = True)
lsc_odr   = odr.ODR(lsc_data, lsc_model, beta0)
lsc_out   = lsc_odr.run()

这些点描述了一个曲率半径为4.5，圆锥常数为-0.8的圆锥曲线。我的代码没有正常工作：通过ODR，代码返回的R值是8（初始值），而不是4.5。有人知道我的代码哪里出问题了吗？

谢谢你的帮助！

Answer 1

你提供的关于Z的数据被忽略了。实际上，你是在计算Z，让它总是符合你定义的隐含方程，无论你传入什么参数。

def f_3(beta, M):

    r = calc_r(M[0],M[1])
    Z = M[2]

    return r**2 - 2*beta[0]*Z + (1+K)*Z**2

...
lsc_data  = odr.Data(row_stack([X, Y, Z]), y=1)

这次运行的结果是R = 4.34911251 +- 0.30341252，看起来和你的预期是吻合的。