matlab做横截面回归,matlab - 将横截面表面轮廓拟合到通用的已知公式以获得系数并对表面进行数学建模 - 堆栈内存溢出...

假设您有一个离散的r数组，并且该数组的每个值z(r) 。 您想拟合曲线以估计非球面透镜的参数。 我将使用lmfit提到这里展现做到这一点使用python的一种方式。

导入用于此的模块：

import numpy as np

import matplotlib.pyplot as plt

from lmfit import Model, Parameters

定义非球面镜的功能：

def asphere_complete(x, r0, k, a2, a4, a6, a8, a10, a12):

r_squared = x ** 2.

z_even_r = r_squared * (a2 + (r_squared * (a4 + r_squared * (a6 + r_squared * (a8 + r_squared * (a10 + (r_squared * a12)))))))

square_root_term = 1 - (1 + k) * ((x / r0) ** 2)

zg = (x ** 2) / (r0 * (1 + np.sqrt(square_root_term)))

return z_even_r + zg

由于您未提供任何数据，因此我将使用以下内容创建一些示例数据，包括人为噪声：

def generate_dummy_data(x, asphere_parameters, noise_sigma, seed=12345):

np.random.seed(seed)

return asphere_complete(x, **asphere_parameters) + noise_sigma * np.random.randn(x.shape[0])

以下函数进行拟合并绘制结果曲线：

def fit_asphere(r, z, fit_parameters):

# create two subplots to plot the original data and the fit in one plot and the residual in another

fig, axarr = plt.subplots(1, 2, figsize=(10, 5))

fit_plot = axarr[0]

residuum_plot = axarr[1]

# configure first plot:

fit_plot.set_xlabel("r")

fit_plot.set_ylabel("z")

fit_plot.grid()

# configure second plot:

residuum_plot.set_xlabel("r")

residuum_plot.set_ylabel("$\Delta$z")

residuum_plot.grid()

# plot original data

fit_plot.plot(r, z, label="Input")

# create an lmfit model and the parameters

function_model = Model(asphere_complete)

# The fitting procedure may throw ValueErrors, if the radicand gets negative

try:

result = function_model.fit(z, fit_parameters, x=r)

# To plot the resulting curve remove the parameters which were just used for the constraints

opt_parameters = dict(result.values)

opt_parameters.pop('r_max', None)

opt_parameters.pop('radicand', None)

# calculate z-values of fitted curve:

z_fitted = asphere_complete(r, **opt_parameters)

# calculate residual values

z_residual = z - z_fitted

# plot fit and residual:

fit_plot.plot(r, z_fitted, label="Fit")

residuum_plot.plot(r, z_residual, label="Residual")

# legends:

fit_plot.legend(loc="best")

residuum_plot.legend(loc="best")

print(result.fit_report())

except ValueError as val_error:

print("Fit Failed: ")

print(val_error)

要设置示例数据的参数，我使用lmfit的Parameters对象：

if __name__ == "__main__":

parameters_dummy = Parameters()

parameters_dummy.add('r0', value=-34.4)

parameters_dummy.add('k', value=-0.98)

parameters_dummy.add('a2', value=0)

parameters_dummy.add('a4', value=-9.67e-9)

parameters_dummy.add('a6', value=1.59e-10)

parameters_dummy.add('a8', value=-5.0e-12)

parameters_dummy.add('a10', value=0)

parameters_dummy.add('a12', value=-1.0e-19)

创建示例数据：

r = np.linspace(0, 35, 1000)

z = generate_dummy_data(r, parameters_dummy, 0.00001)

使用lmfit而不是scipy的curve_fit是，平方根的radicand可能变为负数。 我们需要确保：

为此，我们需要提到定义的约束在这里 。 让我们开始定义我们要在拟合中使用的参数。 基本半径的添加非常简单：

parameters = Parameters()

parameters.add('r0', value=-30, vary=True)

为了服从不等式，请添加变量radicand ，该变量不得小于零。 而不是让k参与拟合常态，而是使r0取决于r0 ， r_max和radicand 。 我们需要使用r_max因为对于最大r ，不等式最成问题。 解决k的不等式导致

用作下面的expr 。 我使用布尔标志来打开/关闭约束：

keep_radicand_safe = True

if keep_radicand_safe:

r_max = np.max(r)

parameters.add('r_max', r_max, vary=False)

parameters.add('radicand', value=0.98, vary=True, min=0)

parameters.add('k', expr='(r0/r_max)**2*(1-radicand)-1')

else:

parameters.add('k', value=-0.98, vary=True)

其余参数可以直接添加：

parameters.add('a2', value=0, vary=False)

parameters.add('a4', value=0, vary=True)

parameters.add('a6', value=0, vary=True)

parameters.add('a8', value=0, vary=True)

parameters.add('a10', value=0, vary=False)

parameters.add('a12', value=0, vary=True)

现在我们准备开始并获得我们的结果：

fit_asphere(r, z, parameters)

plt.show()

在控制台上，您应该看到输出：

[[Variables]]

r0: -34.3999435 +/- 6.1027e-05 (0.00%) (init = -30)

r_max: 35 (fixed)

radicand: 0.71508611 +/- 0.09385813 (13.13%) (init = 0.98)

k: -0.72477176 +/- 0.09066656 (12.51%) == '(r0/r_max)**2*(1-radicand)-1'

a2: 0 (fixed)

a4: 7.7436e-07 +/- 2.7872e-07 (35.99%) (init = 0)

a6: 2.5547e-10 +/- 6.3330e-11 (24.79%) (init = 0)

a8: -4.9832e-12 +/- 1.7115e-14 (0.34%) (init = 0)

a10: 0 (fixed)

a12: -9.8670e-20 +/- 2.0716e-21 (2.10%) (init = 0)

使用我上面使用的数据，如果将keep_radicand_safe设置为False ，则应该看到拟合失败。