Szacowanie przedziałów ufności wokół filtra Kalmana

Question

Pracowałem nad wdrożeniem filtru Kalmana w celu wyszukiwania anomalii w dwuwymiarowym zestawie danych. Bardzo podobny do doskonałego postu, który znalazłem tutaj. W następnym kroku chciałbym przewidzieć przedziały ufności (na przykład 95% pewności dla wartości podłogi i sufitu) dla tego, co przewiduję, że następne wartości będą spadać. Więc oprócz poniższej linii, chciałbym być potrafi wygenerować dwie dodatkowe linie, które zapewniają 95% pewności, że następna wartość będzie powyżej poziomu podłogi lub poniżej pułapu.

Zakładam, że będę chciał użyć macierzy kowariancji niepewności (P), która jest zwracana z każdą prognozą generowaną przez filtr Kalmana, ale nie jestem pewien, czy to prawda. Wszelkie wskazówki lub odniesienia do tego, jak to zrobić, będą mile widziane!

kalman 2d filter in python

Kod na stanowisku powyżej generuje zestaw pomiarów w czasie, stosuje się filtr Kalmana, aby złagodzić efekty.

import numpy as np 
import matplotlib.pyplot as plt 

def kalman_xy(x, P, measurement, R, 
       motion = np.matrix('0. 0. 0. 0.').T, 
       Q = np.matrix(np.eye(4))): 
    """ 
Parameters:  
x: initial state 4-tuple of location and velocity: (x0, x1, x0_dot, x1_dot) 
P: initial uncertainty convariance matrix 
measurement: observed position 
R: measurement noise 
motion: external motion added to state vector x 
Q: motion noise (same shape as P) 
""" 
return kalman(x, P, measurement, R, motion, Q, 
       F = np.matrix(''' 
        1. 0. 1. 0.; 
        0. 1. 0. 1.; 
        0. 0. 1. 0.; 
        0. 0. 0. 1. 
        '''), 
       H = np.matrix(''' 
        1. 0. 0. 0.; 
        0. 1. 0. 0.''')) 

def kalman(x, P, measurement, R, motion, Q, F, H): 
    ''' 
    Parameters: 
    x: initial state 
    P: initial uncertainty convariance matrix 
    measurement: observed position (same shape as H*x) 
    R: measurement noise (same shape as H) 
    motion: external motion added to state vector x 
    Q: motion noise (same shape as P) 
    F: next state function: x_prime = F*x 
    H: measurement function: position = H*x 

    Return: the updated and predicted new values for (x, P) 

    See also http://en.wikipedia.org/wiki/Kalman_filter 

    This version of kalman can be applied to many different situations by 
    appropriately defining F and H 
    ''' 
    # UPDATE x, P based on measurement m  
    # distance between measured and current position-belief 
    y = np.matrix(measurement).T - H * x 
    S = H * P * H.T + R # residual convariance 
    K = P * H.T * S.I # Kalman gain 
    x = x + K*y 
    I = np.matrix(np.eye(F.shape[0])) # identity matrix 
    P = (I - K*H)*P 

    # PREDICT x, P based on motion 
    x = F*x + motion 
    P = F*P*F.T + Q 

    return x, P 

def demo_kalman_xy(): 
    x = np.matrix('0. 0. 0. 0.').T 
    P = np.matrix(np.eye(4))*1000 # initial uncertainty 

    N = 20 
    true_x = np.linspace(0.0, 10.0, N) 
    true_y = true_x**2 
    observed_x = true_x + 0.05*np.random.random(N)*true_x 
    observed_y = true_y + 0.05*np.random.random(N)*true_y 
    plt.plot(observed_x, observed_y, 'ro') 
    result = [] 
    R = 0.01**2 
    for meas in zip(observed_x, observed_y): 
     x, P = kalman_xy(x, P, meas, R) 
     result.append((x[:2]).tolist()) 
    kalman_x, kalman_y = zip(*result) 
    plt.plot(kalman_x, kalman_y, 'g-') 
    plt.show() 

demo_kalman_xy() 

Answer 1

2D uogólnienie 1-sigma interval jest elipsą ufności, który charakteryzuje się wzorem (x-mx).T P^{-1}.(x-mx)==1 z x jest parametrem 2D-wektor, mx oznacza 2D lub ośrodek elipsy i P^{-1} odwrotność macierzy kowariancji. Zobacz, jak narysować jedną z nich: answer. Podobnie jak interwały sigma, powierzchnia elips odpowiada stałemu prawdopodobieństwu, że prawdziwa wartość leży w obrębie. Skalowanie za pomocą współczynnika n (skalowanie długości interwału lub promieni elipsy) pozwala uzyskać wyższą pewność. Zauważ, że czynniki n mają różne prawdopodobieństwa w jednym i dwóch wymiarach:

|`n` | 1D-Intverval | 2D Ellipse | 
================================== 
    1 | 68.27%  | 39.35%  
    2 | 95.5%  | 86.47% 
    3 | 99.73%  | 98.89% 

Obliczenie tych wartości w 2D jest nieco zaangażowane i niestety nie mam publiczne odniesienie do niej.

Answer 2

Jeśli chcesz, aby przedział 95% przewidywał, że następna wartość spadnie, potrzebujesz interwału przewidywania, a nie przedziału ufności (http://en.wikipedia.org/wiki/Prediction_interval).

Dla danych 2-D (3-D), można znaleźć półosi w elipsie (elipsoidalnej) obliczając wartości własne macierzy kowariancji danych i dostosowując rozmiar półosi do rachunku dla niezbędnego prawdopodobieństwa przewidywania.

Wyświetl kod Prediction ellipse and prediction ellipsoid w celu obliczenia elipsy przewidującej 95% lub elipsoidy. To może pomóc w obliczeniu przewidywanej elipsy dla danych.

Answer 3

Ponieważ twoja statystyka pochodzi oczywiście z próbki, prawdopodobieństwo, że statystyczna populacja jest większa od odchylenia standardowego 2 sigma wynosi 0,5. W związku z tym rozważałbym znaczenie rozważenia, czy masz dobrą prognozę wartości, której spodziewasz się, że następny miernik będzie niższy z prawdopodobieństwem 0.95, jeśli nie zastosowałeś wyższego współczynnika ufności 2x odchylenie standardowe. Wielkość tego współczynnika będzie zależeć od wielkości próby użytej do określenia prawdopodobieństwa 0,5 populacji. Im mniejsza jest wielkość próbki użytej do wyprowadzenia macierzy kowariancji, tym większy jest współczynnik do wyliczenia prawdopodobieństwa 0,95, statystyczna populacja 0,95 jest mniejsza od statystycznej statystyki próbki.