我一直在使用Matlab，但我的愿景是最终将所有分析工作转移到Python上进行，因为它是一种真正的编程语言，还有其他几个原因。
我目前尝试解决的问题是对复杂数据进行最小二乘法拟合。作为一名工程师，我们经常处理复杂阻抗，我正在尝试使用曲线拟合将简单电路模型拟合到测量数据上。
阻抗方程如下：
Z(w) = 1/( 1/R + j*w*C ) + j*w*L
然后我尝试找出R、C和L的值，以使得找到最小二乘曲线。
我尝试使用优化包，例如optimize.curve_fit或optimize.leastsq，但它们无法处理复数。
然后我尝试使我的残差函数返回复杂数据的幅度，但这也没有成功。
参考unutbu的答案，在函数残差中取平方和并不必要，因为leastsq无需关心数字是实数还是复数，只需表达为1D数组，并保留功能关系的完整性。
以下是替换后的残差函数：
def residuals(params, w, Z):
R, C, L = params
diff = model(w, R, C, L) - Z
z1d = np.zeros(Z.size*2, dtype = np.float64)
z1d[0:z1d.size:2] = diff.real
z1d[1:z1d.size:2] = diff.imag
return z1d
唯一需要更改的是答案，2013年的答案是：[2.96564781, 1.99929516, 4.00106534]。
相对于unutbu答案的误差降低了一个数量级以上。
最终目标是减少模型和观测值之间平方差的绝对值之和，其中包括Z。
abs(((model(w, *params) - Z)**2).sum())
我的原始回答建议将leastsq应用于一个残差函数，该函数返回一个标量，表示实数和虚数差的平方和：
def residuals(params, w, Z):
R, C, L = params
diff = model(w, R, C, L) - Z
return diff.real**2 + diff.imag**2
Mike Sulzer建议使用返回浮点数向量的剩余函数。
以下是使用这些残差函数的结果比较：
from __future__ import print_function
import random
import numpy as np
import scipy.optimize as optimize
j = 1j
def model1(w, R, C, L):
Z = 1.0/(1.0/R + j*w*C) + j*w*L
return Z
def model2(w, R, C, L):
Z = 1.0/(1.0/R + j*w*C) + j*w*L
# make Z non-contiguous and of a different complex dtype
Z = np.repeat(Z, 2)
Z = Z[::2]
Z = Z.astype(np.complex64)
return Z
def make_data(R, C, L):
N = 10000
w = np.linspace(0.1, 2, N)
Z = model(w, R, C, L) + 0.1*(np.random.random(N) + j*np.random.random(N))
return w, Z
def residuals(params, w, Z):
R, C, L = params
diff = model(w, R, C, L) - Z
return diff.real**2 + diff.imag**2
def MS_residuals(params, w, Z):
"""
https://dev59.com/3nLYa4cB1Zd3GeqPY4cO#20104454 (Mike Sulzer)
"""
R, C, L = params
diff = model(w, R, C, L) - Z
z1d = np.zeros(Z.size*2, dtype=np.float64)
z1d[0:z1d.size:2] = diff.real
z1d[1:z1d.size:2] = diff.imag
return z1d
def alt_residuals(params, w, Z):
R, C, L = params
diff = model(w, R, C, L) - Z
return diff.astype(np.complex128).view(np.float64)
def compare(*funcs):
fmt = '{:15} | {:37} | {:17} | {:6}'
header = fmt.format('name', 'params', 'sum(residuals**2)', 'ncalls')
print('{}\n{}'.format(header, '-'*len(header)))
fmt = '{:15} | {:37} | {:17.2f} | {:6}'
for resfunc in funcs:
# params, cov = optimize.leastsq(resfunc, p_guess, args=(w, Z))
params, cov, infodict, mesg, ier = optimize.leastsq(
resfunc, p_guess, args=(w, Z),
full_output=True)
ssr = abs(((model(w, *params) - Z)**2).sum())
print(fmt.format(resfunc.__name__, params, ssr, infodict['nfev']))
print(end='\n')
R, C, L = 3, 2, 4
p_guess = 1, 1, 1
seed = 2013
model = model1
np.random.seed(seed)
w, Z = make_data(R, C, L)
assert np.allclose(model1(w, R, C, L), model2(w, R, C, L))
print('Using model1')
compare(residuals, MS_residuals, alt_residuals)
model = model2
print('Using model2')
compare(residuals, MS_residuals, alt_residuals)
Using model1
name | params | sum(residuals**2) | ncalls
------------------------------------------------------------------------------------
residuals | [ 2.86950167 1.94245378 4.04362841] | 9.41 | 89
MS_residuals | [ 2.85311972 1.94525477 4.04363883] | 9.26 | 29
alt_residuals | [ 2.85311972 1.94525477 4.04363883] | 9.26 | 29
Using model2
name | params | sum(residuals**2) | ncalls
------------------------------------------------------------------------------------
residuals | [ 2.86590332 1.9326829 4.0450271 ] | 7.81 | 483
MS_residuals | [ 2.85422448 1.94853383 4.04333851] | 9.78 | 754
alt_residuals | [ 2.85422448 1.94853383 4.04333851] | 9.78 | 754
看起来使用哪个剩余函数可能取决于模型函数。考虑到model1和model2的相似性，我无法解释结果的差异。
diff.view('float')更改为diff.astype(np.complex128).view(np.float64)来解决。（我已经编辑了上面的答案以反映这一点。）非连续性问题不会发生，因为diff是两个数组的差异：diff = model(w, R, C, L) - Z。NumPy将始终创建一个包含model(w, R, C, L) - Z结果的连续数组，并且Python将变量名diff绑定到它。 - unutbumodel2的residuals找到的参数比与model1本身配对的任何残差函数都更适合于model1。正如你所指出的，model2返回的值比model1的分辨率低，因此对我来说这是没有意义的，结果应该更好才对。也许这只是运气？尽管如此，我仍然希望知道如何在使用model1时从leastsq中获得相同或更好的结果。 - unutbu