refactorizing code for readability.

Former-commit-id: e62f4ee [formerly e62f4ee [formerly 16ac32e]]
Former-commit-id: 3dc300030e7d4db27f3dc5c7c4453bf43984a9af
Former-commit-id: 1353b1e

Author: compops

python/example1-lgss.py

@@ -1,99 +1,51 @@
-##############################################################################
-#
-# Example of state estimation in a LGSS model 
-# using particle filters and Kalman filters
-#
-# Copyright (C) 2017 Johan Dahlin < liu (at) johandahlin.com.nospam >
-#
-# This program is free software; you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation; either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License along
-# with this program; if not, write to the Free Software Foundation, Inc.,
-# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-#
-##############################################################################
+# State estimation in a LGSS model using particle and Kalman filters
 
 from __future__ import print_function, division
-
-# Import libraries
 import matplotlib.pylab as plt
 import numpy as np
 
-# Import helpers
 from helpers.dataGeneration import generateData
 from helpers.stateEstimation import particleFilter, kalmanFilter
 
 # Set the random seed to replicate results in tutorial
 np.random.seed(10)
 
-
-#=============================================================================
 # Define the model
-#=============================================================================
-
-# Here, we use the following model
-#
-# x[t + 1] = phi * x[t] + sigmav * v[t]
-# y[t] = x[t] + sigmae * e[t]
-#
-# where v[t] ~ N(0, 1) and e[t] ~ N(0, 1)
-
-# Set the parameters of the model (phi, sigmav, sigmae)
-theta = np.zeros(3)
-theta[0] = 0.75
-theta[1] = 1.00
-theta[2] = 0.10
-
-# Set the number of time steps to simulate
-T = 250
-
-# Set the initial state
+# x[t + 1] = phi * x[t] + sigmav * v[t],    v[t] ~ N(0, 1)
+# y[t] = x[t] + sigmae * e[t],              e[t] ~ N(0, 1)
+
+# Set the parameters of the model theta=(phi, sigmav, sigmae), T, x_0
+parameters = np.zeros(3)
+parameters[0] = 0.75
+parameters[1] = 1.00
+parameters[2] = 0.10
+noObservations = 250
 initialState = 0
 
-
-#=============================================================================
 # Generate data
-#=============================================================================
-
-x, y = generateData(theta, T, initialState)
+state, observations = generateData(parameters, noObservations, initialState)
 
-# Plot the measurement
+# Plot data
 plt.subplot(3, 1, 1)
-plt.plot(y, color='#1B9E77', linewidth=1.5)
+plt.plot(observations, color='#1B9E77', linewidth=1.5)
 plt.xlabel("time")
 plt.ylabel("measurement")
 
-# Plot the states
 plt.subplot(3, 1, 2)
-plt.plot(x, color='#D95F02', linewidth=1.5)
+plt.plot(state, color='#D95F02', linewidth=1.5)
 plt.xlabel("time")
 plt.ylabel("latent state")
 
+# State estimation using particle filter with 100 particles
+xHatPF, _ = particleFilter(observations, parameters, 100, initialState)
 
-#=============================================================================
-# State estimation
-#=============================================================================
-
-# Using N = 100 particles and plot the estimate of the latent state
-xHatFilteredParticleFilter, _ = particleFilter(y, theta, 100, initialState)
-
-# Using the Kalman filter
-xHatFilteredKalmanFilter = kalmanFilter(y, theta, initialState, 0.01)
+# State estimation using the Kalman filter
+xHatKF = kalmanFilter(observations, parameters, initialState, 0.01)
 
+# Plot state estimate
 plt.subplot(3, 1, 3)
-plt.plot(xHatFilteredKalmanFilter[1:T] - xHatFilteredParticleFilter[0:T-1], color='#7570B3', linewidth=1.5)
+plt.plot(xHatKF[1:noObservations] - xHatPF[0:noObservations-1], color='#7570B3', linewidth=1.5)
 plt.xlabel("time")
 plt.ylabel("difference in estimate")
-plt.show()
 
-##############################################################################
-# End of file
-##############################################################################
+plt.show()

python/example2-lgss.py

@@ -1,130 +1,72 @@
-##############################################################################
-#
-# Example of particle Metropolis-Hastings in a LGSS model.
-#
-# Copyright (C) 2017 Johan Dahlin < liu (at) johandahlin.com.nospam >
-#
-# This program is free software; you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation; either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License along
-# with this program; if not, write to the Free Software Foundation, Inc.,
-# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-#
-##############################################################################
+# Parameter estimation using particle Metropolis-Hastings in a LGSS model.
 
 from __future__ import print_function, division
-
-# Import libraries
 import matplotlib.pylab as plt
 import numpy as np
 
-# Import helpers
 from helpers.dataGeneration import generateData
 from helpers.stateEstimation import particleFilter, kalmanFilter
 from helpers.parameterEstimation import particleMetropolisHastings
 
-# Set the random seed
+# Set the random seed to replicate results in tutorial
 np.random.seed(10)
 
-
-#=============================================================================
 # Define the model
-#=============================================================================
-
-# Here, we use the following model
-#
-# x[t + 1] = phi * x[t] + sigmav * v[t]
-# y[t] = x[t] + sigmae * e[t]
-#
-# where v[t] ~ N(0, 1) and e[t] ~ N(0, 1)
-
-# Set the parameters of the model (phi, sigmav, sigmae)
-theta = np.zeros(3)
-theta[0] = 0.75
-theta[1] = 1.00
-theta[2] = 0.10
-
-# Set the number of time steps to simulate
-T = 250
-
-# Set the initial state
+# x[t + 1] = phi * x[t] + sigmav * v[t],    v[t] ~ N(0, 1)
+# y[t] = x[t] + sigmae * e[t],              e[t] ~ N(0, 1)
+
+# Set the parameters of the model theta = (phi, sigmav, sigmae), T, x_0
+parameters = np.zeros(3)
+parameters[0] = 0.75
+parameters[1] = 1.00
+parameters[2] = 0.10
+noObservations = 250
 initialState = 0
 
-
-#=============================================================================
-# Generate data
-#=============================================================================
-
-x, y = generateData(theta, T, initialState)
-
-
-#=============================================================================
-# Parameter estimation using PMH
-#=============================================================================
-
-# The inital guess of the parameter
+# Settings for PMH
 initialPhi = 0.50
-
-# No. particles in the particle filter ( choose noParticles ~ T )
-noParticles = 500
-
-# The length of the burn-in and the no. iterations of PMH 
-# ( noBurnInIterations < noIterations )
+noParticles = 500           # Use noParticles ~ noObservations
 noBurnInIterations = 1000
 noIterations = 5000
-
-# The standard deviation in the random walk proposal
 stepSize = 0.10
 
-# Run the PMH algorithm
-res = particleMetropolisHastings(y, initialPhi, theta, noParticles, initialState,
-                                 particleFilter, noIterations, stepSize)
+# Generate data
+state, observations = generateData(parameters, noObservations, initialState)
 
+# Run the PMH algorithm
+phiTrace = particleMetropolisHastings(
+    observations, initialPhi, parameters, noParticles, 
+    initialState, particleFilter, noIterations, stepSize)
 
-#=============================================================================
 # Plot the results
-#=============================================================================
-
 noBins = int(np.floor(np.sqrt(noIterations - noBurnInIterations)))
 grid = np.arange(noBurnInIterations, noIterations, 1)
-resPhi = res[noBurnInIterations:noIterations]
+phiTrace = phiTrace[noBurnInIterations:noIterations]
 
-# Plot the parameter posterior estimate
-# Solid black line indicate posterior mean
+# Plot the parameter posterior estimate (solid black line = posterior mean)
 plt.subplot(3, 1, 1)
-plt.hist(resPhi, noBins, normed=1, facecolor='#7570B3')
+plt.hist(phiTrace, noBins, normed=1, facecolor='#7570B3')
 plt.xlabel("phi")
 plt.ylabel("posterior density estimate")
-plt.axvline(np.mean(resPhi), color='k')
+plt.axvline(np.mean(phiTrace), color='k')
 
-# Plot the trace of the Markov chain after burn-in
-# Solid black line indicate posterior mean
+# Plot the trace of the Markov chain after burn-in (solid black line = posterior mean)
 plt.subplot(3, 1, 2)
-plt.plot(grid, resPhi, color='#7570B3')
+plt.plot(grid, phiTrace, color='#7570B3')
 plt.xlabel("iteration")
 plt.ylabel("phi")
-plt.axhline(np.mean(resPhi), color='k')
+plt.axhline(np.mean(phiTrace), color='k')
 
 # Plot the autocorrelation function
 plt.subplot(3, 1, 3)
-macf = np.correlate(resPhi - np.mean(resPhi), resPhi - np.mean(resPhi), mode='full')
-macf = macf[macf.size/2:]
+macf = np.correlate(phiTrace - np.mean(phiTrace), phiTrace - np.mean(phiTrace), mode='full')
+idx = int(macf.size/2)
+macf = macf[idx:]
 macf = macf[0:100]
 macf /= macf[0]
 grid = range(len(macf))
 plt.plot(grid, macf, color='#7570B3')
 plt.xlabel("lag")
 plt.ylabel("ACF of phi")
 
-
-##############################################################################
-# End of file
-##############################################################################
+plt.show()