非线性拟合正弦曲线失败答案

【问题标题】：Failure of non linear fit to sine curve非线性拟合正弦曲线失败
【发布时间】：2013-08-13 17:14:43
【问题描述】：

给定一些生成的二维玩具数据，我一直在尝试拟合正弦曲线的幅度、频率和相位。（代码在最后）

为了获得三个参数的估计值，我首先执行 FFT。我使用 FFT 中的值作为实际频率和相位的初始猜测，然后对其进行拟合（逐行）。我编写了我的代码，以便我输入我希望频率在哪个 FFT 箱中，这样我就可以检查拟合是否正常工作。但是有一些非常奇怪的行为。如果我的输入 bin 是 3.1（一个非整数 bin，所以 FFT 不会给我正确的频率），那么拟合效果很好。但是如果输入 bin 是 3（因此 FFT 输出准确的频率），那么我的拟合失败了，我正试图了解原因。

这是我将输入箱（在 X 和 Y 方向）分别设为 3.0 和 2.1 时的输出：

（右边的图是数据-拟合）

这是我将输入箱设为 3.0 和 2.0 时的输出：

问题：为什么我输入曲线的准确频率时非线性拟合会失败？

代码：

#! /usr/bin/python

# For the purposes of this code, it's easier to think of the X-Y axes as transposed, 
# so the X axis is vertical and the Y axis is horizontal

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as optimize
import itertools
import sys

PI = np.pi

# Function which accepts paramters to define a sin curve
# Used for the non linear fit    
def sineFit(t, a, f, p):
   return a * np.sin(2.0 * PI * f*t + p)

xSize    = 18
ySize    = 60
npt      = xSize * ySize

# Get frequency bin from user input
xFreq    = float(sys.argv[1])
yFreq    = float(sys.argv[2])

xPeriod  = xSize/xFreq
yPeriod  = ySize/yFreq

# arrays should be defined here

# Generate the 2D sine curve
for jj in range (0, xSize):
   for ii in range(0, ySize):
      sineGen[jj, ii] = np.cos(2.0*PI*(ii/xPeriod + jj/yPeriod))

# Compute 2dim FFT as well as freq bins along each axis
fftData  = np.fft.fft2(sineGen)
fftMean  = np.mean(fftData)
fftRMS   = np.std(fftData)
xFreqArr = np.fft.fftfreq(fftData.shape[1]) # Frequency bins along x
yFreqArr = np.fft.fftfreq(fftData.shape[0]) # Frequency bins along y

# Find peak of FFT, and position of peak
maxVal = np.amax(np.abs(fftData))
maxPos = np.where(np.abs(fftData) == maxVal)

# Iterate through peaks in the FFT 
# For this example, number of loops will always be only one

prevPhase = -1000
for col, row in itertools.izip(maxPos[0], maxPos[1]):

   # Initial guesses for fit parameters from FFT
   init_phase  = np.angle(fftData[col,row])
   init_amp    = 2.0 * maxVal/npt
   init_freqY  = yFreqArr[col]
   init_freqX  = xFreqArr[row]

   cntr  = 0
   if prevPhase == -1000:
      prevPhase = init_phase

   guess = [init_amp, init_freqX, prevPhase]
   # Fit each row of the 2D sine curve independently
   for rr in sineGen:   
      (amp, freq, phs), pcov = optimize.curve_fit(sineFit, xDat, rr, guess)
      # xDat is an linspace array, containing a list of numbers from 0 to xSize-1

      # Subtract fit from original data and plot
      fitData     = sineFit(xDat, amp, freq, phs)
      sub1        = rr - fitData

      # Plot
      fig1 = plt.figure()
      ax1  = fig1.add_subplot(121)
      p1,  = ax1.plot(rr, 'g')
      p2,  = ax1.plot(fitData, 'b')
      plt.legend([p1,p2], ["data", "fit"])

      ax2  = fig1.add_subplot(122)
      p3,  = ax2.plot(sub1)
      plt.legend([p3], ['residual1'])

      fig1.tight_layout()

      plt.show()
      cntr += 1
      prevPhase = phs # Update guess for phase of sine curve

【问题讨论】：

我不可能阅读这么多代码（没有 cmets！）并找到错误...我建议将您的程序分解成碎片，一个构建所有这些数组，另一个构建合适的，另一个制作情节的...然后分别测试它们，找出错误所在并重新发布，以便我们可以查看它并使用一小段代码。
您已经发布了代码，但这不是一个最小的工作示例。即使有 cmets，为什么没有定义 xDat 或 sineGen？
@Hooked - 对不起......为了让我的代码更小，我只是用 cmets 代替了定义。我会重新编辑它！

标签： python scipy curve-fitting

【解决方案1】：

我已尝试将您问题的重要部分提炼成这个答案。

首先，尝试拟合单个数据块，而不是数组。一旦你确信你的模型足够了，你就可以继续前进。
您的合身度只会与您的模型一样好，如果您继续使用非“正弦”的东西，则需要相应地进行调整。
拟合是一门“艺术”，因为初始条件可以极大地改变误差函数的收敛性。此外，您的拟合可能不止一个最小值，因此您经常需要担心您提出的解决方案的独特性。

虽然您的 FFT 想法走在了正确的轨道上，但我认为您的实现并不完全正确。下面的代码应该是一个很棒的玩具系统。它生成f(x) = a0*sin(a1*x+a2) 类型的随机数据。有时随机的初始猜测会起作用，有时会失败。但是，使用 FFT 猜测频率收敛应该始终适用于该系统。示例输出：

import numpy as np
import pylab as plt
import scipy.optimize as optimize

# This is your target function
def sineFit(t, (a, f, p)):
    return a * np.sin(2.0*np.pi*f*t + p)

# This is our "error" function
def err_func(p0, X, Y, target_function):
    err = ((Y - target_function(X, p0))**2).sum()
    return err


# Try out different parameters, sometimes the random guess works
# sometimes it fails. The FFT solution should always work for this problem
inital_args = np.random.random(3)

X = np.linspace(0, 10, 1000)
Y = sineFit(X, inital_args)

# Use a random inital guess
inital_guess = np.random.random(3)

# Fit
sol = optimize.fmin(err_func, inital_guess, args=(X,Y,sineFit))

# Plot the fit
Y2 = sineFit(X, sol)
plt.figure(figsize=(15,10))
plt.subplot(211)
plt.title("Random Inital Guess: Final Parameters: %s"%sol)
plt.plot(X,Y)
plt.plot(X,Y2,'r',alpha=.5,lw=10)

# Use an improved "fft" guess for the frequency
# this will be the max in k-space
timestep = X[1]-X[0]
guess_k = np.argmax( np.fft.rfft(Y) )
guess_f = np.fft.fftfreq(X.size, timestep)[guess_k]
inital_guess[1] = guess_f 

# Guess the amplitiude by taking the max of the absolute values
inital_guess[0] = np.abs(Y).max()

sol = optimize.fmin(err_func, inital_guess, args=(X,Y,sineFit))
Y2 = sineFit(X, sol)

plt.subplot(212)
plt.title("FFT Guess          : Final Parameters: %s"%sol)
plt.plot(X,Y)
plt.plot(X,Y2,'r',alpha=.5,lw=10)
plt.show()

【讨论】：

我意识到这是一个“健谈”的评论，但你的数字很漂亮。当我将来显示模型预测与实际数据时，我将不得不开始调整 alpha 设置。谢谢你。
@Engineero 感谢您的补充！在生产中，我使用matplotlibrc 文件来保存我的所有自定义设置（强烈推荐用于图形一致性）。此外，我为标签和绘图打开了 LaTeX 渲染，这让它更有专业的感觉。

【解决方案2】：

问题是由于相位的初始猜测错误，而不是频率。在循环通过 genSine（内循环）的行时，您使用上一行的拟合结果作为下一行的初始猜测，这并不总是有效。如果您从当前行的 fft 确定相位并将其用作初始猜测，则拟合将成功。您可以按如下方式更改内部循环：

for n,rr in enumerate(sineGen):   
    fftx = np.fft.fft(rr)
    fftx = fftx[:len(fftx)/2]
    idx = np.argmax(np.abs(fftx))
    init_phase = np.angle(fftx[idx])
    print fftx[idx], init_phase
    ...

你也需要改变

def sineFit(t, a, f, p):
   return a * np.sin(2.0 * np.pi * f*t + p)

到

def sineFit(t, a, f, p):
   return a * np.cos(2.0 * np.pi * f*t + p)

因为 phase=0 意味着 fft 的虚部为零，因此该函数类似于余弦。

顺便说一句。您上面的示例仍然缺少 sineGen 和 xDat 的定义。

【讨论】：

【解决方案3】：

根据http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html，在不了解您的大部分代码的情况下：

(amp2, freq2, phs2), pcov = optimize.curve_fit(sineFit, tDat, 
                                                     sub1, guess2)

应该变成：

(amp2, freq2, phs2), pcov = optimize.curve_fit(sineFit, tDat, 
                                                         sub1, p0=guess2)

假设 tDat 和 sub1 是 x 和 y，这应该可以解决问题。但是，再一次，理解如此复杂的代码是相当困难的，它有如此多的相互关联的变量，而根本没有 cmets。代码应始终自下而上构建，这意味着当单个代码不工作时，您不会进行拟合循环，在代码工作以适应非嘈杂示例之前，您不会添加噪声......好好运！

【讨论】：

我已经对其进行了编辑以添加 cmets... 抱歉发布这样的代码，如果还有什么需要更改的请告诉我！
感谢cmets，我稍后会尝试看看。同时，您可以尝试查看没有噪音或任何花哨的单一拟合是否适用于您的代码的拟合部分。
这几乎就是我所做的一切......以前的代码很混乱，但我在我的问题中发布的情节正是这种情况。没有噪音，并且是一排二维正弦曲线，如果我提供的频率恰好在一个 FFT 箱中，则拟合似乎失败了。
查看下面的新评论。

【解决方案4】：

“没什么花哨的”我的意思是删除与合身无关的所有内容，并做一个简化的模拟示例，例如：

import numpy as np
import scipy.optimize as optimize

def sineFit(t, a, f, p):
       return a * np.sin(2.0 * np.pi * f*t + p)


# Create array of x and y with given parameters
x = np.asarray(range(100))
y = sineFit(x, 1, 0.05, 0)

# Give a guess and fit, printing result of the fitted values
guess = [1., 0.05, 0.]
print optimize.curve_fit(sineFit, x, y, guess)[0]

这个结果就是答案：

[1.    0.05   0.]

但如果你改变猜测不要太多，就够了：

# Give a guess and fit, printing result of the fitted values
guess = [1., 0.06, 0.]
print optimize.curve_fit(sineFit, x, y, guess)[0]

结果给出了荒谬的错误数字：

[ 0.00823701  0.06391323 -1.20382787]

你能解释一下这种行为吗？

【讨论】：

这是答案还是问题？
这是试图缩小帖子问题的范围。我试图说明这可能（或可能不是）是问题的一部分。

【解决方案5】：

您可以将curve_fit 与一系列三角函数一起使用，这些函数通常非常健壮并且可以通过增加项的数量来调整到您需要的精度...这是一个示例：

from scipy import sin, cos, linspace
def f(x, a0,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,
            c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12):
    return a0 + s1*sin(1*x) +  c1*cos(1*x) \
              + s2*sin(2*x) +  c2*cos(2*x) \
              + s3*sin(3*x) +  c3*cos(3*x) \
              + s4*sin(4*x) +  c4*cos(4*x) \
              + s5*sin(5*x) +  c5*cos(5*x) \
              + s6*sin(6*x) +  c6*cos(6*x) \
              + s7*sin(7*x) +  c7*cos(7*x) \
              + s8*sin(8*x) +  c8*cos(8*x) \
              + s9*sin(9*x) +  c9*cos(9*x) \
             + s10*sin(9*x) + c10*cos(9*x) \
             + s11*sin(9*x) + c11*cos(9*x) \
             + s12*sin(9*x) + c12*cos(9*x)

from scipy.optimize import curve_fit
pi/2. / (x.max() - x.min())
x_norm *= norm_factor
popt, pcov = curve_fit(f, x_norm, y)
x_fit = linspace(x_norm.min(), x_norm.max(), 1000)
y_fit = f(x_fit, *popt)
plt.plot( x_fit/x_norm, y_fit )

【讨论】：