Stima intervalli di confidenza attorno al filtro Kalman

Question

Ho lavorato per implementare un filtro Kalman per cercare anomalie in un insieme di dati bidimensionale. Molto simile all'eccellente post che ho trovato qui. Come passo successivo, mi piacerebbe prevedere gli intervalli di confidenza (ad esempio il 95% di confidenza per i valori del pavimento e del soffitto) per quello che prevedo che i prossimi valori cadranno. Quindi oltre alla linea qui sotto, mi piacerebbe essere in grado di generare due linee aggiuntive che rappresentano una confidenza del 95% che il valore successivo sarà sopra il pavimento o sotto il soffitto.

Suppongo di voler utilizzare la matrice di covarianza di incertezza (P) restituita con ogni previsione generata dal filtro di Kalman ma non sono sicuro che sia corretta. Qualunque guida o riferimento su come farlo sarebbe molto apprezzato!

kalman 2d filter in python

Il codice nel post genera sopra di una serie di misure nel tempo e utilizza un filtro di Kalman per lisciare i risultati.

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

La generalizzazione 2D del 1-sigma interval è l'ellisse fiducia che è caratterizzato dall'equazione (x-mx).T P^{-1}.(x-mx)==1, con x essendo il parametro 2D-vettore, mx 2D significa o centro dell'ellisse e P^{-1} la matrice di covarianza inversa. Vedi questo answer su come disegnarne uno. Come gli intervalli sigma, l'area delle ellissi corrisponde a una probabilità fissa che il vero valore sia all'interno. Ridimensionando il fattore n (ridimensionando la lunghezza dell'intervallo oi raggi dell'ellisse) è possibile raggiungere una maggiore sicurezza. Si noti che i fattori n hanno diverse probabilità in una e due dimensioni:

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

Il calcolo di questi valori in 2D è un po 'coinvolto e purtroppo non ho un riferimento pubblico ad esso.

Answer 2

Se si desidera un intervallo del 95% per prevedere i valori successivi, allora si desidera un intervallo di previsione e non un intervallo di confidenza (http://en.wikipedia.org/wiki/Prediction_interval).

Per i dati 2-D (3-D), i semi-assi dell'ellisse (ellissoide) possono essere trovati calcolando gli autovalori della matrice di covarianza dei dati e regolando la dimensione dei semi-assi sul conto per la probabilità di previsione necessaria.

Vedere Prediction ellipse and prediction ellipsoid per un codice Python per calcolare l'ellisse o ellissoide di previsione del 95%. Questo potrebbe aiutare a calcolare l'ellisse di previsione per i dati.

Answer 3

Poiché la statistica è ovviamente derivata da un campione, la probabilità che la statistica della popolazione sia maggiore della deviazione standard 2 sigma è 0,5. Pertanto, vorrei contemplare l'importanza di considerare se hai una buona previsione di un valore che ti aspetti che la prossima misura sia inferiore alla probabilità 0,95 se non hai applicato un fattore di confidenza superiore della deviazione standard 2x. L'entità di tale fattore dipenderà dalla dimensione del campione utilizzata per ricavare la probabilità di 0.5 popolazione. Più piccola è la dimensione del campione utilizzata per ricavare la matrice di covarianza, più grande è il fattore per ricavare la probabilità di 0,95 che la popolazione 0,95 di statistica è inferiore alla statistica campionaria su fattorizzata.