rlabbe / filterpy

Python Kalman filtering and optimal estimation library. Implements Kalman filter, particle filter, Extended Kalman filter, Unscented Kalman filter, g-h (alpha-beta), least squares, H Infinity, smoothers, and more. Has companion book 'Kalman and Bayesian Filters in Python'.
MIT License
3.31k stars 617 forks source link

Corrected shape for self.sigmas_h #168

Closed Prokhozhijj closed 5 years ago

Prokhozhijj commented 6 years ago

This correction is related to the bug #129

rlabbe commented 6 years ago

This is not correct. You never provided source code for #129 to duplicate that bug.

With this change, none of the unit tests for the UKF pass. The unit test test, among other things, that the behavior of the UKF is identical to the KF when the problem is linear, so this is very strong evidence that the code is substantially correct.

Now, it may be the code doesn't handle some edge case that you are using, or your problem formulation is wrong. I need to see source code to duplicate your problem.

Ideally, any submission should include a description of what is failing, along with source code to duplicate it, and, best of all, a unit test that fails with the current implementation and that would pass with a correct implementation. Without that I am left guessing as to what problem you may be having. I want my code to be both correct and robust across any problem formulation, but I cannot do this without the necessary information.

Prokhozhijj commented 6 years ago

Of course you are right. I should be more attentive to preliminary code testing. I am sorry for caused inconveniences. Here is my code that leads to the issue #129.

from filterpy.kalman import UnscentedKalmanFilter as UKF
from filterpy.kalman import (unscented_transform, MerweScaledSigmaPoints,
                             JulierSigmaPoints, SimplexSigmaPoints)
from filterpy.common import Q_discrete_white_noise

import pandas as pd
import numpy as np

def read_data():
    return np.random.rand(1000)

def get_filter(dim_x, dt, Q_mult, R, z0):
    def f_Fv(x, dt):
        """ Process transition model with speed only
        """
        A = np.array ([[1, dt],
                       [0, 1]])
        return A.dot(x)

    def f_Fa(x, dt):
        """ Process transition model with speed and acceleration
        """
        A = np.array ([[1, dt, dt*dt/2],
                       [0, 1, dt],
                       [0, 0, 1]])
        return A.dot(x)

    def f_Fj(x, dt):
        """ Process transition model with speed, acceleration and jerk
        """
        dt2 = dt * dt
        dt3 = dt2 * dt
        A = np.array ([[1, dt, dt2/2, dt3/6],
                       [0,  1,    dt, dt2/2],
                       [0,  0,     1,    dt],
                       [0,  0,     0,     1]])
        return A.dot(x)

    def f_Hx(x):
        """ Measurement function. Converts state vector x into a measurement
            vector of shape (dim_z = 1).
        """
        return x[0]

    _dim_x = dim_x
    _points = MerweScaledSigmaPoints(n = _dim_x, alpha = 0.001, beta = 2., kappa = 0.)
    _P_mult = 1 
    _z0 = z0
    _Q_mult = Q_mult
    _R = R
    _dt = dt

    # Depending on how much dimensions do we have we need to build process model 
    Q = np.eye(_dim_x)
    x = np.zeros(_dim_x)
    x[0] = _z0

    functions = {2:f_Fv, 3:f_Fa, 4:f_Fj}
    f_F = functions[_dim_x]

    filter = UKF(dim_x = _dim_x, dim_z = 1, dt = _dt, hx = f_Hx, fx = f_F, points = _points)

    filter.x = x

    # process uncertainty (covariance)
    # the more uncertainty, the less we believe in process model
    filter.Q = Q_discrete_white_noise(dim = _dim_x, dt = _dt, var = 0.01)
    filter.Q *= _Q_mult

    # noise measurement (covariance)
    # the less the value, the more we believe in measurement (the price from stock exchange)
    filter.R = _R        

    # uncertainty covariance of starting position
    filter.P *= _P_mult

    return filter

z = read_data()
n = len(z)
start_dt = 150
step_dt = 100
model_noise = 10
process_noise = 10
dim_x = 3

ofilter = get_filter(dim_x = dim_x, dt = 1/start_dt, Q_mult = model_noise, R = process_noise, z0 = z[0])
Xs, _ = ofilter.batch_filter(z)
results = Xs[:, 0]
Prokhozhijj commented 6 years ago

And here is code output:

Traceback (most recent call last): File "ukf_test.py", line 88, in <module> Xs, _ = ofilter.batch_filter(z) File "g:\proj\python\filterpy\filterpy\kalman\UKF.py", line 605, in batch_filter self.update(z, r, UT=UT) File "g:\proj\python\filterpy\filterpy\kalman\UKF.py", line 450, in update zp, self.S = UT(self.sigmas_h, self.Wm, self.Wc, R, self.z_mean, self.residual_z) File "g:\proj\python\filterpy\filterpy\kalman\unscented_transform.py", line 104, in unscented_transform x = np.dot(Wm, sigmas) # dot = \Sigma^n_1 (W[k]*Xi[k]) ValueError: shapes (7,) and (1,7) not aligned: 7 (dim 0) != 1 (dim 0)

rlabbe commented 5 years ago

Your computation for the measurement function should read

def f_Hx(x):
    """ Measurement function. Converts state vector x into a measurement
        vector of shape (dim_z = 1).
    """
    return [x[0]]

As the comment states, it should return a vector. Your version returned x[0] which is a scalar, not a vector. With this code change the filter runs and the output of plt.plot(results) looks reasonable to me, given that you are passing in random data.

Another possible problem is your assignment to R. filter.R = _R is incorrect, since _R is a scalar. filter.R *= _R assigns it the value np.array([[10]]), which is what you want.

Prokhozhijj commented 5 years ago

Thanks! Now it works as expected.

AlexanderDamen commented 4 years ago

Your computation for the measurement function should read

def f_Hx(x):
    """ Measurement function. Converts state vector x into a measurement
        vector of shape (dim_z = 1).
    """
    return [x[0]]

As the comment states, it should return a vector. Your version returned x[0] which is a scalar, not a vector. With this code change the filter runs and the output of plt.plot(results) looks reasonable to me, given that you are passing in random data.

Another possible problem is your assignment to R. filter.R = _R is incorrect, since _R is a scalar. filter.R *= _R assigns it the value np.array([[10]]), which is what you want.

Hey rlabbe, first of all thanks for writing this book. It is really usefull. I had the exact same problem and it took me more than an hour to find this solution. Is it possible to add a check to see if it are vectors or not? I think this would save people a lot of time.