直線データのノイズを除去したい

Question

前提・実現したいこと

y軸の正方向のノイズを含む直線データからノイズを取り除き、直線を取り出すプログラムをつくっています。ノイズの大きさは1%程度です。直線に対して高速フーリエ変換（FFT）は不適なので、移動平均とスプライン補間による平滑化を試みました。しかし、以下のコードでは取り出したい直線をうまく抽出できず困っています。ノイズを含むy=10-x+noise（グラフのlinear, 青）からノイズを含まないy=10-x（グラフのoriginal, 黒）を取り出すのが目標です。改善点、アドバイス、別のアプローチなどがあれば教えていただけると幸いです。
。

以下のサイトを参考にしました。
https://snova301.hatenablog.com/entry/2018/10/07/135233

実行結果

拡大図

該当のソースコード

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d


def make_func(in_x):
    np.random.seed(0)
    out_y = 10 - in_x + 0.3*(np.random.rand(in_x.size)) # 10-x に正方向のみの誤差を加える
    return out_y


def spline_interp(in_x, in_y):
    out_x = np.linspace(np.min(in_x), np.max(in_x), np.size(in_x)*100) # もとのxの個数より多いxを用意
    func_spline = interp1d(in_x, in_y, kind='cubic') # 3次のスプライン曲線
    out_y = func_spline(out_x) # func_splineはscipyオリジナルの型

    return out_x, out_y


def moving_avg(in_x, in_y):
    np_y_conv = np.convolve(in_y, np.ones(3)/float(3), mode='valid')
    out_x_dat = np.linspace(np.min(in_x), np.max(in_x), np.size(np_y_conv))

    return out_x_dat, np_y_conv


def main():
    x1 = np.linspace(0, 10, 100)
    y1 = make_func(x1)

    x2, y2 = spline_interp(x1, y1)

    x3, y3 = moving_avg(x1, y1)
    x4, y4 = spline_interp(x3, y3)
    
    plt.plot(x1, 10-x1, color='k', label='original', alpha=0.7)
    plt.plot(x1, y1, color='b', label='linear', alpha=0.7)
    plt.plot(x2, y2, color='r', label='spline', alpha=0.7)
    plt.plot(x4, y4, color='g', label='average + spline', alpha=0.7)
    plt.legend()
    plt.show()


if __name__ == '__main__':
    main()

補足情報（FW/ツールのバージョンなど）

python 3.8.8

Comments

[jbpb0]

> linear（青）からoriginal（黒）を取り出す

どれ？

[cream_puff]

わかりづらくて申し訳ありません。編集後の内容で伝わりますでしょうか？

[jbpb0]

質問者さんが、グラフの画像をアップロードするのを忘れているのだと思って、質問したのですが

もしかしたら、自分でコードを実行して、表示されるグラフを見ろ、ってこと？

[jbpb0]

ノイズの素性は既知なのでしょうか？
どんな分布の形状(ガウスノイズ？など)かとか、ノイズの大きさはどれくらいかとか

[cream_puff]

申し訳ありませんでした。ノイズの大きさは１％程度、形状は不明です。測定データなのでノイズの大きさはガウス分布に従うと思いますが、ノイズはy軸正方向に現れる可能性が高いです。

[jbpb0]

ノイズの素性が明らかでは無いのなら、データからノイズの分布形状や大きさを推定して、ノイズの影響によるy+方向の平均オフセットを算出して、最小二乗法等の方法で決めた直線を算出したオスセット分y-方向にずらす、みたいなことをやらないといけないような気がします

> ノイズはy軸正方向に現れる可能性が高い

「可能性が高い」ということは、頻度は低いけれどy-方向のノイズも有り得る、ということですか？

もしそうなら、質問を編集して、そのように書いてください
ppaulさんは、「ノイズが絶対に正の値である」ことを前提に回答を書いてます
マイナスのノイズも(頻度は低くても)有り得るなら、話が変わります

[ozwk]

ノイズが正方向なのではなくて、測定データに正のオフセットが乗っているのでは?

Answer 1

ノイズが真の値に与える影響をどのように仮定するかによりますが、以下のように仮定します。

y ~ uniform(min = intercept + coef * x , max = intercept + coef * x + shift)

intercept + coef * x の部分は真の値を表し、一様分布の最小値とします。ノイズがshiftであり、何らかの固定値とします。実際の観測値は、真の値とノイズを加えた値の間のどこかに一様分布に従って出現するものとします。

上記のようなモデルにおいてintercept,coef,shiftを推計することができれば求める真値を得ることができます。方法としてはstanを使ってはどうでしょうか。実装は以下になります。

#### 学習結果を安定させるため、1)係数の正負を反転、2)y1の最小値がゼロから始まるように調整

y1_ = -y1 + np.max(y1)

import pystan  #### pythonからStanをキックするためのモジュール

#### モデル本体の定義

stan_code = """
data { 
    int N  ; // data数
    real X[N] ; // 説明変数
    real Y[N] ; // 被説明変数
    }
    
parameters {
    real coef ;           // 係数
    real intercept ;      // 切片
    real<lower=0> shift ; // 正方向へのシフト値
}

model {
    for (i in 1:N) {
        Y[i] ~ uniform(coef * X[i] + intercept ,coef * X[i] + intercept + shift) ;
    }
}

"""

model = pystan.StanModel(model_code=stan_code)  ### モデルのコンパイル

#### stan投入用データの編集
stan_data = {'N' : len(x1),
             'X' : x1,
             'Y' : y1_}

#### モデルの学習
result = model.sampling(data=stan_data,iter=6000,chains=1,seed=1234,warmup=3000)

#### 学習結果からcoef,intercept,shiftを取得して元の値に戻す
Coef = np.mean(result.extract()['coef']) * -1
Intercept = np.mean(result.extract()['intercept']) + np.max(y1)
Shift = np.mean(result.extract()['shift'])

#### 推計結果の算出
predict = Coef*x1 + Intercept - Shift

上記は、あくまでもノイズが一様分布の場合に限ります。他の分布によるのであれば、stan部分を変更する必要がありますし、場合によってはpreproccessing部分は必要ありません。

Answer 2

移動平均の横軸の値の算出方法がおかしいので、修正例を書いときます
ついでに、平均の長さを簡単に変えられるようにしました

def moving_avg(in_x, in_y, ave_len):
    if ave_len % 2 == 0: 
        ave_len = ave_len + 1
    np_y_conv = np.convolve(in_y, np.ones(ave_len)/float(ave_len), mode='valid')
    #out_x_dat = np.linspace(np.min(in_x), np.max(in_x), np.size(np_y_conv))
    ave_len_h = int((ave_len - 1) / 2)
    out_x_dat = np.linspace(in_x[ave_len_h], in_x[-ave_len_h-1], np.size(np_y_conv))
    return out_x_dat, np_y_conv


def main():
    x1 = np.linspace(0, 10, 100)
    y1 = make_func(x1)
    x2, y2 = spline_interp(x1, y1)
    x3, y3 = moving_avg(x1, y1, 21)
    x4, y4 = spline_interp(x3, y3)
    plt.plot(x1, 10-x1, color='k', label='original', alpha=0.7)
    plt.plot(x1, y1, color='b', label='linear', alpha=0.7)
    #plt.plot(x2, y2, color='r', label='spline', alpha=0.7)
    plt.plot(x3, y3, color='r', label='average', alpha=0.7)
    #plt.plot(x4, y4, color='g', label='average + spline', alpha=0.7)
    plt.legend()
    plt.show()

ただし、あくまでも、移動平均のみの修正で、ノイズによるオフセットの補正は入ってません

平均の長さを長くすれば、そこそこ直線っぽくはなります
ただし、直線と分かっているのなら、最小二乗法で直線にフィットする方がいいです
移動平均を勧めているわけではありません

Comments

[cream_puff]

ありがとうございます。参考にさせていただきます。

Answer 3

ノイズが絶対に正の値であるということであれば、以下の方法があります。

import numpy as np
import matplotlib.pyplot as plt

def make_func(in_x, noise):
    np.random.seed(0)
    out_y = 10 - in_x + noise*(np.random.rand(in_x.size)) # 10-x に正方向のみの誤差を加える
    return out_y

def main(number, noise):
    x1 = np.linspace(0, 10, number+1)
    y1 = make_func(x1, noise)
    xy1 = np.polyfit(x1,y1,1)
    y2 = np.poly1d(xy1)(x1)
    y3 = y2 + min(y1-y2)
    plt.plot(x1, 10-x1, color='k', label='original', alpha=0.7)
    plt.plot(x1, y1, color='r', label='noise', alpha=0.7)
    plt.plot(x1, y2, color='b', label='linear', alpha=0.7)
    plt.plot(x1, y3, color='g', label='adjust', alpha=0.7)
    plt.legend()
    plt.show()

number=100、noise=2.0だと少しずれます。

number=1000、noise=2.0というようにサンプル数を増やすと、ほぼ元の関数に一致します。

Comments

[cream_puff]

ありがとうございます。参考にさせていただきます。

Answer 4

生方向のみのノイズが乗っている、とするなら、なにをしたところでもとの信号の復元は無理です
普通は正負両方向のノイズが乗っているという前提で移動平均を掛けてノイズを除去するものです

Answer 5

ノイズの量や性質にもよりますが、こんな感じで直線を近似するのはどうでしょうか？

    xy1 = np.polyfit(x1,y1,1)
    y5 = np.poly1d(xy1)(x1)

    plt.plot(x1, y5, color='r', label='fit', alpha=0.7)

Comments

[cream_puff]

早速ありがとうございます。ノイズを含んだまま直線でフィッティングすると取り出したい直線をy軸正方向に平行移動した直線 y = 10 - x + (ノイズの平均値) が出力されます。ノイズを取り除いた y = 10 - x が出力するようにしたいと考えております。

[ppaul]

cream_puffさんは中心極限定理をご存知ですか。
ご存じでなければ
https://ja.wikipedia.org/wiki/%E4%B8%AD%E5%BF%83%E6%A5%B5%E9%99%90%E5%AE%9A%E7%90%86
をお読み下さい。

[cream_puff]

ppaulさんありがとうございます。読ませていただきました。標本平均はガウス分布に従うということは、フィッティング結果は y = 10 - x + (ガウス分布) となるということでしょうか。