5

私はここで私が理解していないことが何であるかを理解しようとしています。

私はhttp://www.scipy.org/Cookbook/FittingDataをフォローしており、正弦波に適合させようとしています。本当の問題は、回転する宇宙船で素晴らしい正弦波を作る衛星磁力計のデータです。データセットを作成し、それを適合させて入力を復元しようとしています。

これが私のコードです：

import numpy as np
from scipy import optimize

from scipy.optimize import curve_fit, leastsq

import matplotlib.pyplot as plt


class Parameter:
    def __init__(self, value):
            self.value = value

    def set(self, value):
            self.value = value

    def __call__(self):
            return self.value

def fit(function, parameters, y, x = None):
    def f(params):
        i = 0
        for p in parameters:
            p.set(params[i])
            i += 1
        return y - function(x)

    if x is None: x = np.arange(y.shape[0])
    p = [param() for param in parameters]
    return optimize.leastsq(f, p, full_output=True, ftol=1e-6, xtol=1e-6)

# generate a perfect data set (my real data have tiny error)
def mysine(x, a1, a2, a3):
    return a1 * np.sin(a2 * x + a3)

xReal = np.arange(500)/10.
a1 = 200.
a2 = 2*np.pi/10.5  # omega, 10.5 is the period
a3 = np.deg2rad(10.) # 10 degree phase offset
yReal = mysine(xReal, a1, a2, a3)

# plot the real data
plt.figure(figsize=(15,5))
plt.plot(xReal, yReal, 'r', label='Real Values')

# giving initial parameters
amplitude = Parameter(175.)
frequency = Parameter(2*np.pi/8.)
phase = Parameter(0.0)

# define your function:
def f(x): return amplitude() * np.sin(frequency() * x + phase())

# fit! (given that data is an array with the data to fit)
out = fit(f, [amplitude, frequency, phase], yReal, xReal)
period = 2*np.pi/frequency()
print amplitude(), period, np.rad2deg(phase())

xx = np.linspace(0, np.max(xReal), 50)
plt.plot( xx, f(xx) , label='fit')
plt.legend(shadow=True, fancybox=True)

これはこのプロットを作ります： ここに画像の説明を入力してください

の復元された適合パラメータは、[44.2434221897 8.094832581 -61.6204033699]私が始めたものとは似ていません。

私が理解していないことや間違っていることについて何か考えはありますか？

scipy.__version__
'0.10.1'

編集：1つのパラメータを修正することが提案されました。上記の例では、振幅を固定してnp.histogram(yReal)[1][-1]も許容できない出力が生成されます。適合：[175.0 8.31681375217 6.0]別の適合方法を試す必要がありますか？どの提案？

ここに画像の説明を入力してください

4

2 に答える 2

2

From what I can see from playing a bit with leastsq (without fancy stuff from the cookbook, just plain direct calls to leastsq --- and by the way, full_output=True is your friend here), is that it's very hard to fit all three of the amplitude, frequency and phase in one go. On the other hand, if I fix the amplitude and fit the frequency and phase, it works; if I fix the frequency and fit the amplitude and phase, it works too.

There is more than one way out here. What might be the simplest one --- if you are sure you only have one sine wave (and this is easy to check with the Fourier transform), then you know the frequency from just the distance between consecutive maxima of your signal. Then fit the two remaining parameters.

If what you have is a mixture of several harmonics, well, again, Fourier transform will tell you that.

于 2012-11-15T21:19:19.410 に答える
2

これは、Zhenyaのアイデアのいくつかを実装するコードです。それは使用しています

yhat = fftpack.rfft(yReal)
idx = (yhat**2).argmax()
freqs = fftpack.rfftfreq(N, d = (xReal[1]-xReal[0])/(2*pi))
frequency = freqs[idx]

データの主な頻度を推測し、

amplitude = yReal.max()

振幅を推測します。

import numpy as np
import scipy.optimize as optimize
import scipy.fftpack as fftpack
import matplotlib.pyplot as plt
pi = np.pi
plt.figure(figsize = (15, 5))

# generate a perfect data set (my real data have tiny error)
def mysine(x, a1, a2, a3):
    return a1 * np.sin(a2 * x + a3)

N = 5000
xmax = 10
xReal = np.linspace(0, xmax, N)
a1 = 200.
a2 = 2*pi/10.5  # omega, 10.5 is the period
a3 = np.deg2rad(10.) # 10 degree phase offset
print(a1, a2, a3)
yReal = mysine(xReal, a1, a2, a3) + 0.2*np.random.normal(size=len(xReal))

yhat = fftpack.rfft(yReal)
idx = (yhat**2).argmax()
freqs = fftpack.rfftfreq(N, d = (xReal[1]-xReal[0])/(2*pi))
frequency = freqs[idx]

amplitude = yReal.max()
guess = [amplitude, frequency, 0.]
print(guess)
(amplitude, frequency, phase), pcov = optimize.curve_fit(
    mysine, xReal, yReal, guess)

period = 2*pi/frequency
print(amplitude, frequency, phase)

xx = xReal
yy = mysine(xx, amplitude, frequency, phase)
# plot the real data
plt.plot(xReal, yReal, 'r', label = 'Real Values')
plt.plot(xx, yy , label = 'fit')
plt.legend(shadow = True, fancybox = True)
plt.show()

収量

(200.0, 0.5983986006837702, 0.17453292519943295)   # (a1, a2, a3)
[199.61981404516041, 0.61575216010359946, 0.0]     # guess
(200.06145097308041, 0.59841420869261097, 0.17487141943703263) # fitted parameters

fftを使用することにより、頻度の推測はすでに最終的に適合したパラメーターにかなり近いことに注意してください。

パラメータを修正する必要はないようです。頻度の推測を実際の値に近づけることによりoptimize.curve_fit、合理的な答えに収束することができます。

于 2012-11-15T22:13:38.567 に答える