Extended Kalman Filter Non-Linear Regression Implementation R, Correct?

Question

I am working on the Kalman Filter and its applications. I tried to implement a model for a nonlinear regression problem of the form:

$$ y = \exp(-(X\beta)) + q,\quad q \sim \mathcal{N}(0, I) $$

using the Extended Kalman-Filter. Please note that I assume that the derivation of the equations to be known, as their explanation would go beyond the scope of this question. Further model specifications can be found in the code. My implementation does not seem to be totally wrong, but for large values of $y$, the prediction is a bit out of control. Any suggestions about potential for improvement, or validation of correctness are greatly appreciated.

Here is my implementation in R:

# number of samples to generate
n <- 500

# regressors with intercept
X <- cbind(rep(1, n), rnorm(n))

# true state
beta <- c(0.1, 1)

# observations wit normal error
y <- exp(-(X %*% beta)) + rnorm(n, sd=1)

# plot data, non linear!
plot( 1 + X[,2], y)

# state transition matrix is the identity:
stateTransition <- diag(c(1,1))

# 2 states
n_state <- ncol(stateTransition)    

# init variable to store the states in loop
state <- matrix(0, ncol=n_state, nrow=n_measure)

# process noise is the identity
processNoise <- diag(c(1,1))

# measurement noise is 1
measurementNoise <- 1
S <- diag(rep(1, n_state))

# iterate
for(i in 2:n){
    # predict
    m_ <- stateTransition %*% state[i-1, ]
    S_ <- stateTransition %*% S %*% t(stateTransition) + processNoise

    # update equations, drop = false preseves the structure of the matrix, otherwise R would
    # turn it into a vector
    v <- y[i] - exp(-(X[i,,drop=FALSE] %*% m_))

    # Jakobian of the model
    H <- -X[i,,drop=FALSE] * exp(-(X[i,,drop=FALSE] %*% m_))[1]

    # computations according to model
    S <- H  %*% S_ %*% t(H) + measurementNoise
    K <- S_ %*% t(H) %*% solve(S)
    m <- m_ + K %*% v
    S <- S_ - K %*% S %*% t(K)
    # save state
    state[i,] <- m
}

# generate predictions according to estimated state
pred <- numeric(n)
for(i in 1:n){
    pred[i] <- exp(-(X[i,]%*%state[i,]))
}

# plot results
plot(y)
lines(pred, col='red')

As we can see, for moderately high values of $y$ the prediction is stable; however for the outliers, the prediction explodes.

Answer 1

1. In your model, beta is fixed, so you should set processNoise <- diag(c(0, 0)). By specifying non-zero variances for the state, you'd assume that the state (beta) is stochastic.
2. As your initial guess of beta is arbitrary, you should set the initial uncertainty of beta to a relatively big value, or use an exact initialisation. In this specific example, S <- diag(rep(1, n_state)) seems to be sufficient, as it is equal to the variance of the measurements.
3. The line S <- S_ - K %*% S %*% t(K) seems to be wrong, it should be S <- S_ - K %*% H %*% S_
4. You refer to n_measure which is not specified, should be n in order for the code to run properly.
5. You're not using the first observation. To fix this, specify m <- 0 just before the for loop, let the loop start at i = 1, and replace state[i-1, ] by m.
6. You're talking about predictions, but strictly speaking, these are filtered signals.