2017-04-12 20 views
14

[編集] @Claudioの答えは、私に外れ値を除外する方法についての本当に良いヒントを与えます。私は自分のデータでカルマンフィルタを使い始めたいと思っています。そこで以下のサンプルデータを変更して極端ではない微妙なバリエーションノイズが発生するようにしました(私も同様です)。他の誰かが私のデータにPyKalmanをどのように使うべきかについての指示を私に与えることができれば、それは素晴らしいものになるでしょう。 場所データのPythonでカルマンフィルタを使用するには?

ロボットプロジェクトの場合、私はカメラで空気中の凧を追跡しようとしています。私はPythonでプログラミングしています。以下のようなノイズの多い場所の結果を貼り付けました(すべてのアイテムにはdatetimeオブジェクトが含まれていますが、わかりやすくするためにそれらを残しました)。

[   # X  Y 
    {'loc': (399, 293)}, 
    {'loc': (403, 299)}, 
    {'loc': (409, 308)}, 
    {'loc': (416, 315)}, 
    {'loc': (418, 318)}, 
    {'loc': (420, 323)}, 
    {'loc': (429, 326)}, # <== Noise in X 
    {'loc': (423, 328)}, 
    {'loc': (429, 334)}, 
    {'loc': (431, 337)}, 
    {'loc': (433, 342)}, 
    {'loc': (434, 352)}, # <== Noise in Y 
    {'loc': (434, 349)}, 
    {'loc': (433, 350)}, 
    {'loc': (431, 350)}, 
    {'loc': (430, 349)}, 
    {'loc': (428, 347)}, 
    {'loc': (427, 345)}, 
    {'loc': (425, 341)}, 
    {'loc': (429, 338)}, # <== Noise in X 
    {'loc': (431, 328)}, # <== Noise in X 
    {'loc': (410, 313)}, 
    {'loc': (406, 306)}, 
    {'loc': (402, 299)}, 
    {'loc': (397, 291)}, 
    {'loc': (391, 294)}, # <== Noise in Y 
    {'loc': (376, 270)}, 
    {'loc': (372, 272)}, 
    {'loc': (351, 248)}, 
    {'loc': (336, 244)}, 
    {'loc': (327, 236)}, 
    {'loc': (307, 220)} 
] 

まず、異常値を手作業で計算してから、リアルタイムでそれらを単にデータから削除することを考えました。次に、私はカルマンフィルタと、それが具体的にどのようにノイズの多いデータを滑らかにすることになっているかについて読んでいます。 検索の結果、PyKalman libraryが見つかりました。私はカルマンフィルタの用語全体で少し失われていたので、私はwikiとカルマンフィルタのいくつかの他のページを読みました。私はカルマンフィルタの一般的な考え方を得ていますが、私は自分のコードにどのように適用すればいいか分かりません。 PyKalman docs

私は、次の例を見つけました:

>>> from pykalman import KalmanFilter 
>>> import numpy as np 
>>> kf = KalmanFilter(transition_matrices = [[1, 1], [0, 1]], observation_matrices = [[0.1, 0.5], [-0.3, 0.0]]) 
>>> measurements = np.asarray([[1,0], [0,0], [0,1]]) # 3 observations 
>>> kf = kf.em(measurements, n_iter=5) 
>>> (filtered_state_means, filtered_state_covariances) = kf.filter(measurements) 
>>> (smoothed_state_means, smoothed_state_covariances) = kf.smooth(measurements) 

次のように私は単に私自身の観察のための観測を置換:

from pykalman import KalmanFilter 
import numpy as np 
kf = KalmanFilter(transition_matrices = [[1, 1], [0, 1]], observation_matrices = [[0.1, 0.5], [-0.3, 0.0]]) 
measurements = np.asarray([(399,293),(403,299),(409,308),(416,315),(418,318),(420,323),(429,326),(423,328),(429,334),(431,337),(433,342),(434,352),(434,349),(433,350),(431,350),(430,349),(428,347),(427,345),(425,341),(429,338),(431,328),(410,313),(406,306),(402,299),(397,291),(391,294),(376,270),(372,272),(351,248),(336,244),(327,236),(307,220)]) 
kf = kf.em(measurements, n_iter=5) 
(filtered_state_means, filtered_state_covariances) = kf.filter(measurements) 
(smoothed_state_means, smoothed_state_covariances) = kf.smooth(measurements) 

をそれは私にどんな意味のあるデータを与えるものではありません。

>>> smoothed_state_means 
array([[-235.47463353, 36.95271449], 
     [-354.8712597 , 27.70011485], 
     [-402.19985301, 21.75847069], 
     [-423.24073418, 17.54604304], 
     [-433.96622233, 14.36072376], 
     [-443.05275258, 11.94368163], 
     [-446.89521434, 9.97960296], 
     [-456.19359012, 8.54765215], 
     [-465.79317394, 7.6133633 ], 
     [-474.84869079, 7.10419182], 
     [-487.66174033, 7.1211321 ], 
     [-504.6528746 , 7.81715451], 
     [-506.76051587, 8.68135952], 
     [-510.13247696, 9.7280697 ], 
     [-512.39637431, 10.9610031 ], 
     [-511.94189431, 12.32378146], 
     [-509.32990832, 13.77980587], 
     [-504.39389762, 15.29418648], 
     [-495.15439769, 16.762472 ], 
     [-480.31085928, 18.02633612], 
     [-456.80082586, 18.80355017], 
     [-437.35977492, 19.24869224], 
     [-420.7706184 , 19.52147918], 
     [-405.59500937, 19.70357845], 
     [-392.62770281, 19.8936389 ], 
     [-388.8656724 , 20.44525168], 
     [-361.95411607, 20.57651509], 
     [-352.32671579, 20.84174084], 
     [-327.46028214, 20.77224385], 
     [-319.75994982, 20.9443245 ], 
     [-306.69948771, 21.24618955], 
     [-287.03222693, 21.43135098]]) 

は私よりも明るく魂が正しい方向に私にいくつかのヒントや例を与えることができる:たとえば、smoothed_state_meansは、次のようになり?すべてのヒントは大歓迎です!

+1

フィルタが必要な場合もありますが、カルマンフィルタが必要かどうかわかりません。あなたがカルマンフィルタを必要としていると確信していない限り、私はここで使用するフィルタリングの種類について尋ねることをお勧めします:http://dsp.stackexchange.com/ –

+0

あなたの質問に答えません。 3シグマ以外の値を削除すると、投稿されたノイズのある値はすべて削除されます。 – Ben

+0

私の(微妙な)理解では、カルマンフィルタは、(不完全な)物理的/数学的モデルと実際の(ノイズの多い)測定値との間の不一致を調整します。 - あなたの問題の声明では、ポジションの予測モデルを認識できないので、カルマンフィルタがあなたに役立つかどうか疑問に思います。 – gboffi

答えて

12

TL; DR、下のコードと画像を参照してください。

私はカルマンフィルタがあなたのアプリケーションではかなりうまくいくと思っていますが、カイトのダイナミクス/物理学についてもう少し考えが必要です。

this webpageを強くお勧めします。私は著者との関係や知識はありませんが、カルマンフィルタのフィルターを丸ごと鳴らそうと約1日を費やしました。

短く; (すなわち、状態と入力を知っていれば、将来の状態を予測することができます)、それが本当の状態であると推定するためにシステムについて知っているものを組み合わせる最適な方法を提供します。測定ノイズ」の対象となっている

  1. 測定(:(あなたがそれを記述したページに表示すべての行列代数によって考慮される)巧妙なビットはそれが最適にあなたが持っている2つの情報を組み合わせた方法です(つまり、センサーが完璧ではない)

  2. ダイナミクス(状態が「プロセスノイズ」の影響を受けやすい入力を受けてどのように進化していると考えているか) 。

あなたはこれらの各(共同差異行列RQを経由して、それぞれ)上にあるかどうかはわかり指定し、カルマンゲインはあなたのモデル(すなわちを信じるべきかどの位決定あなたの現在のあなたの国の推定値)、そしてあなたの測定値をどれくらい信じるべきか。

あなたのカイトのシンプルなモデルを構築しましょう。以下で私が提案するのは、非常にシンプルなモデルです。カイトのダイナミクスについてもっと知っている方がいいかもしれません。

便宜のために、我々は状態ベクトルで書くことができる四つの状態を有し、(その3次元の方向性を持っている、明らかに簡素化、本物の凧は、拡張体である)のは、粒子として凧を扱うみましょう:

X = [X、x_dot、Y、y_dot]、xおよびyは位置であり、_dotするのは、これらの方向の各々における速度である

。あなたの質問から、私たちは測定ベクトルで書くことができる2つの(潜在的に騒々しい)の測定は、そこにあると仮定しています:

Z = [X、Y]は、

を私たちが書くダウンすることができ、測定を行列(Hhereを論じ、そしてpykalmanライブラリでobservation_matrices):

Z = H X =>H = [1、0、0、0]、 [0、0、1、0]]

次に、システムのダイナミクスについて説明する必要があります。ここで私は、外力は作用せず、カイトの動きにはダンピングがないと仮定します(知識があれば、より良くできるかもしれませんが、これは外力を効果的に扱い、未知の/モデル化されていない外乱として効果的に扱います)。

この場合、以前のサンプル「K-1」の状態の関数として現在のサンプル「K」に我々の状態のそれぞれのためのダイナミクスは以下のように与えられる:

X(K)=×( K-1)+ DT * x_dot(K-1)

x_dot(K)= x_dot(K-1)

Y(K)= Y(K-1)+ DT * y_dot(K- 1)

y_dot(K)= y_dot(K-1) "DT" はティムで

e-step。現在の位置と速度に基づいて(x、y)位置が更新され、速度は変更されないままであると仮定する。ユニットが与えられていないとすれば、上記の式、すなわちposition_units/sample_intervalの単位で "dt"を省略できるような速度単位があると言うことができます(私はあなたの測定サンプルが一定の間隔であると仮定します)。我々のようにダイナミクスマトリックスへのこれらの4つの方程式をまとめることができる(Fここで説明し、transition_matricesライブラリpykalmanで): Fxの

X(K)= (K-1)=>F = [[1,1,0,0]、[0,1,0,0]、[0,0,1,1]、[0、0、0、1]]である。

ここで、PythonでKalmanフィルタを使用することができます。コードから変更:

(青赤、x位置であるy座標であり、x軸は単にサンプル数である)、それはノイズを拒否する合理的仕事をして示す以下を作製
from pykalman import KalmanFilter 
import numpy as np 
import matplotlib.pyplot as plt 
import time 

measurements = np.asarray([(399,293),(403,299),(409,308),(416,315),(418,318),(420,323),(429,326),(423,328),(429,334),(431,337),(433,342),(434,352),(434,349),(433,350),(431,350),(430,349),(428,347),(427,345),(425,341),(429,338),(431,328),(410,313),(406,306),(402,299),(397,291),(391,294),(376,270),(372,272),(351,248),(336,244),(327,236),(307,220)]) 

initial_state_mean = [measurements[0, 0], 
         0, 
         measurements[0, 1], 
         0] 

transition_matrix = [[1, 1, 0, 0], 
        [0, 1, 0, 0], 
        [0, 0, 1, 1], 
        [0, 0, 0, 1]] 

observation_matrix = [[1, 0, 0, 0], 
         [0, 0, 1, 0]] 

kf1 = KalmanFilter(transition_matrices = transition_matrix, 
        observation_matrices = observation_matrix, 
        initial_state_mean = initial_state_mean) 

kf1 = kf1.em(measurements, n_iter=5) 
(smoothed_state_means, smoothed_state_covariances) = kf1.smooth(measurements) 

plt.figure(1) 
times = range(measurements.shape[0]) 
plt.plot(times, measurements[:, 0], 'bo', 
     times, measurements[:, 1], 'ro', 
     times, smoothed_state_means[:, 0], 'b--', 
     times, smoothed_state_means[:, 2], 'r--',) 
plt.show() 

enter image description here

あなたは上記のプロットを見て、それがあまりにもでこぼこに見えると思うと仮定します。どうやってそれを修正できますか?カルマンフィルタ上述したように二つの情報に作用している:

  • システムダイナミクス(及び状態の現在の推定値)(私達の状態の2つのxおよびyのこの場合)

    1. 測定

    上記のモデルで取得されたダイナミクスは非常に簡単です。文字通り、彼らは現在の速度(現実的で妥当な方法で)によって位置が更新され、速度は一定であると言います(これは明らかに物理的に真ではありませんが、速度がゆっくりと変化するという直感を捉えています)。

    推定状態がより滑らかでなければならないと考える場合、これを達成する1つの方法は、私たちのダイナミクスよりも測定値に対する信頼度が低い(すなわち、我々のstate_covarianceに対してより高いobservation_covariance)。示されるようにem_vars設定が観測共分散の再推定を回避するために必要とされる上記のコードの端から開始

    、前回推定observation covariance 10Xの値を修正し、(here参照)

    生成
    kf2 = KalmanFilter(transition_matrices = transition_matrix, 
            observation_matrices = observation_matrix, 
            initial_state_mean = initial_state_mean, 
            observation_covariance = 10*kf1.observation_covariance, 
            em_vars=['transition_covariance', 'initial_state_covariance']) 
    
    kf2 = kf2.em(measurements, n_iter=5) 
    (smoothed_state_means, smoothed_state_covariances) = kf2.smooth(measurements) 
    
    plt.figure(2) 
    times = range(measurements.shape[0]) 
    plt.plot(times, measurements[:, 0], 'bo', 
         times, measurements[:, 1], 'ro', 
         times, smoothed_state_means[:, 0], 'b--', 
         times, smoothed_state_means[:, 2], 'r--',) 
    plt.show() 
    

    以下のプロット(点としての測定値、点線の状態推定値)。違いはやや微妙ですが、うまくいけばスムーズだと分かります。

    enter image description here

    あなたがオンラインでこのフィットフィルターを使用したい場合は最後に、あなたはfilter_updateメソッドを使用して行うことができます。 smoothメソッドはバッチ測定にしか適用できないため、これはsmoothメソッドではなくfilterメソッドを使用することに注意してください。 hereもっと:以下

    time_before = time.time() 
    n_real_time = 3 
    
    kf3 = KalmanFilter(transition_matrices = transition_matrix, 
            observation_matrices = observation_matrix, 
            initial_state_mean = initial_state_mean, 
            observation_covariance = 10*kf1.observation_covariance, 
            em_vars=['transition_covariance', 'initial_state_covariance']) 
    
    kf3 = kf3.em(measurements[:-n_real_time, :], n_iter=5) 
    (filtered_state_means, filtered_state_covariances) = kf3.filter(measurements[:-n_real_time,:]) 
    
    print("Time to build and train kf3: %s seconds" % (time.time() - time_before)) 
    
    x_now = filtered_state_means[-1, :] 
    P_now = filtered_state_covariances[-1, :] 
    x_new = np.zeros((n_real_time, filtered_state_means.shape[1])) 
    i = 0 
    
    for measurement in measurements[-n_real_time:, :]: 
        time_before = time.time() 
        (x_now, P_now) = kf3.filter_update(filtered_state_mean = x_now, 
                 filtered_state_covariance = P_now, 
                 observation = measurement) 
        print("Time to update kf3: %s seconds" % (time.time() - time_before)) 
        x_new[i, :] = x_now 
        i = i + 1 
    
    plt.figure(3) 
    old_times = range(measurements.shape[0] - n_real_time) 
    new_times = range(measurements.shape[0]-n_real_time, measurements.shape[0]) 
    plt.plot(times, measurements[:, 0], 'bo', 
         times, measurements[:, 1], 'ro', 
         old_times, filtered_state_means[:, 0], 'b--', 
         old_times, filtered_state_means[:, 2], 'r--', 
         new_times, x_new[:, 0], 'b-', 
         new_times, x_new[:, 2], 'r-') 
    
    plt.show() 
    

    プロットfilter_update方法を用いて見出さ3点を含むフィルタ方式の性能を示しています。点は測定値であり、破線はフィルタトレーニング期間の状態推定値であり、実線は「オンライン」期間の状態推定値である。

    enter image description here

    そして、(私のラップトップ上の)タイミング情報。

    Time to build and train kf3: 0.0677888393402 seconds 
    Time to update kf3: 0.00038480758667 seconds 
    Time to update kf3: 0.000465154647827 seconds 
    Time to update kf3: 0.000463008880615 seconds 
    
  • +0

    これを外れ値の検出/除去手法と比較することができます。カイトモデルがダイナミクスを仮定していない場合、カルマンフィルタは純粋に測定のノイズがゼロ平均とノーマル分布であると仮定して(平均位置の)最尤推定であると思います。しかし、これは、カイトが急速に下方に移動している場合、次の間隔ではそれ以上ダウンする可能性があるという情報を失うことになります。上で実装された下側では、大きな偏差が完全に無視されるわけではありません(そのためには、「ゲーティングゲーティング」を参照してください)。 – kabdulla

    +0

    精巧な答えをありがとう。それは私に数時間かかりましたが、私は今、カルマンフィルタの考え方がより良くなっていると思います。もう一つの質問。私がスムージングをもっと強くして、ラインが今よりも滑らかになるようにしたいのであれば、どうすればいいでしょうか? – kramer65

    +0

    心配はいりません。あまりに精巧ではないことを願っています。私はカルマンフィルターにはかなり新しいので、簡単に保つようにしていましたが、今私は彼らに熱心ですので、持ち去られたかもしれません!あなたが示唆しているのはかなり単純明快でなければなりません。私はこれについてもう少し詳しく解説します。 – kabdulla

    2

    私がカルマンフィルタを使用して見ることのできるものは、あなたの場合には適切なツールではないかもしれません。

    どうすればいいですか? :

     0  1 
    0 346 226 
    1 346 211 
    2 347 196 
    3 347 180 
    4 350 2165 
    5 355 154 
    6 359 138 
    7 368 120 
    8 374 -830 
    9 346 90 
    10 349 75 
    11 1420 67 
    12 357 64 
    13 358 62 
         0 1 
    0 346 226 
    1 346 211 
    2 347 196 
    3 347 180 
    5 355 154 
    6 359 138 
    7 368 120 
    9 346 90 
    10 349 75 
    12 357 64 
    13 358 62 
    

    は私がから上記のコードを持っているいくつかのより多くのとソースのhereを参照してください:ここでは、出力

    lstInputData = [ 
        [346, 226 ], 
        [346, 211 ], 
        [347, 196 ], 
        [347, 180 ], 
        [350, 2165], ## noise 
        [355, 154 ], 
        [359, 138 ], 
        [368, 120 ], 
        [374, -830], ## noise 
        [346, 90 ], 
        [349, 75 ], 
        [1420, 67 ], ## noise 
        [357, 64 ], 
        [358, 62 ] 
    ] 
    
    import pandas as pd 
    import numpy as np 
    df = pd.DataFrame(lstInputData) 
    print(df) 
    from scipy import stats 
    print (df[(np.abs(stats.zscore(df)) < 1).all(axis=1)]) 
    

    +0

    これは本当に良いアイデアです!私は間違いなくこのトリックを使用するつもりです。これに加えて私はカルマンフィルタを使いたいと思っています。私はカルマンフィルタをどのように使用できるかの例を私にも教えていただけますか? – kramer65

    +2

    @ kramer65私は、カルマンフィルタのフィルタリングを使用する主題があまりにも広すぎてここで議論することはできないと思います。加えて、私はカルマンフィルターの専門家ではないので、あなたが私の答えと一緒に生きていけずそれを受け入れることができないなら、あなたは他の答えを待たなければなりません。あなたが書いたので、私はここで答えました。「すべてのヒントは大歓迎です! :) – Claudio

    関連する問題