fixed conflicts.

Former-commit-id: e39a9fc

Author: Johan Dahlin

README.md

@@ -2,18 +2,22 @@
 
 This code was downloaded from https://github.com/compops/pmh-tutorial and contains the code used to produce the results in the tutorial:
 
-J. Dahlin and T. B. Schön, **Getting started with particle Metropolis-Hastings for inference in nonlinear models**. Pre-print, arXiv:1511:01707, 2017. 
+J. Dahlin and T. B. Schön, **Getting started with particle Metropolis-Hastings for inference in nonlinear models**. Pre-print, arXiv:1511:01707, 2017.
 
-The tutorial is available as a preprint from http://arxiv.org/pdf/1511.01707. An R package is also provided on CRAN with the implementation of the turorial in R. The source code (almost identitical to the code in the subdirectory R/) is found at https://github.com/compops/pmh-tutorial-rpkg.
+The tutorial is available as a preprint from http://arxiv.org/pdf/1511.01707. An R package is also provided on CRAN with the implementation of the tutorial in R. The source code (almost identical to the code in the subdirectory R/) is found at https://github.com/compops/pmh-tutorial-rpkg.
 
 Included material
 --------------
-**r/** This is the main implementation. The complete R code developed and implemented in the tutorial. This code was used to make all the numerical illustrations in the tutorial including the figures and tables. The workspaces for these runs are also provided as a [zip-file in the latest release of the code](https://github.com/compops/pmh-tutorial/releases/latest) to reproduce all the figures in the tutorial. 
+**r/** This is the main implementation. The complete R code developed and implemented in the tutorial. This code was used to make all the numerical illustrations in the tutorial including the figures and tables. The workspaces for these runs are also provided as a [zip-file in the latest release of the code](https://github.com/compops/pmh-tutorial/releases/latest) to reproduce all the figures in the tutorial.
 
 **python/** Code for Python to implement the basic algorithms covered in the tutorial. Implementations for the advanced topics are not provided. Only simple plotting is implemented and no figures or saved data from runs are provided.
 
 **matlab/** Code for MATLAB to implement the basic algorithms covered in the tutorial. Implementations for the advanced topics are not provided. Only simple plotting is implemented and no figures or saved data from runs are provided.
 
+Generalisations
+--------------
+There is source code available for Python that implements some of the generalisations discussed in the tutorial. See the README file under *python/* for more information.
+
 Copyright information
 --------------
 See *LICENSE* for more information.

matlab/README.md

@@ -1,6 +1,6 @@
 # MATLAB code for PMH tutorial
 
-This MATLAB code implements the Kalman filter (KF), particle filter (PF) and particle Metropolis-Hastings (PMH) algorithm for two different dynamical models: a linear Gaussian state-space (LGSS) model and a stochastic volatilty (SV) model. Note that the Kalman filter can only be employed for the first of these two models. The details of the code is described in the tutorial paper available at http://arxiv.org/pdf/1511.01707
+This MATLAB code implements the Kalman filter (KF), particle filter (PF) and particle Metropolis-Hastings (PMH) algorithm for two different dynamical models: a linear Gaussian state-space (LGSS) model and a stochastic volatility (SV) model. Note that the Kalman filter can only be employed for the first of these two models. The details of the code is described in the tutorial paper available at http://arxiv.org/pdf/1511.01707
 
 Note that the MATLAB code in this folder covers the basic implementations in the paper. The notation of the variables has been changed sligthly compared with the tutorial paper to improve readability of the code. However, it should be easy to translate between the two. See the R code in r/ for all the implementations and to recreate the results in the tutorial.
 
@@ -49,4 +49,4 @@ Copyright information
 # with this program; if not, write to the Free Software Foundation, Inc.,
 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 ##############################################################################
-```
+```

python/README.md

@@ -1,6 +1,6 @@
 # Python code for PMH tutorial
 
-This Python code implements the Kalman filter (KF), particle filter (PF) and particle Metropolis-Hastings (PMH) algorithm for two different dynamical models: a linear Gaussian state-space (LGSS) model and a stochastic volatilty (SV) model. Note that the Kalman filter can only be employed for the first of these two models. The details of the code is described in the tutorial paper available at http://arxiv.org/pdf/1511.01707
+This Python code implements the Kalman filter (KF), particle filter (PF) and particle Metropolis-Hastings (PMH) algorithm for two different dynamical models: a linear Gaussian state-space (LGSS) model and a stochastic volatility (SV) model. Note that the Kalman filter can only be employed for the first of these two models. The details of the code is described in the tutorial paper available at http://arxiv.org/pdf/1511.01707
 
 Note that the Python code in this folder covers the basic implementations in the paper. The notation of the variables has been changed sligthly compared with the tutorial paper to improve readability of the code. However, it should be easy to translate between the two. See the R code in r/ for all the implementations and to recreate the results in the tutorial.
 
@@ -26,7 +26,7 @@ Supporting files (helpers/)
 --------------
 **dataGeneration.py** Generates data from a LGSS model.
 
-**parameterEstimation.py** Implements the PMH algorithm for the LGSS model (particleMetropolisHastings) and the SV model (particleMetropolisHastingsSVModel). 
+**parameterEstimation.py** Implements the PMH algorithm for the LGSS model (particleMetropolisHastings) and the SV model (particleMetropolisHastingsSVModel).
 
 **stateEstimation.py** Implements the faPF for the LGSS model (particleFilter), the Kalman filter for the LGSS model (kalmanFilter) and the bPF for the SV model (paticleFilterSVmodel).
 
@@ -35,6 +35,10 @@ Adapting the code for another model
 --------------
 See the discussion in *README.MD* in the directory *r/*.
 
+Generalisations
+--------------
+Some generalisations and improvements of this code is discussed in the tutorial, see the last paragraph in Section 7. Python code for PMH1 and PMH2 is available in the repo [pmh-stco2015](https://github.com/compops/pmh-stco2015), Python code for qPMH2 is availabe in the repo [https://github.com/compops/qpmh2-sysid2015](qpmh2-sysid2015) and Python code for correlated pseudo-marginal Metropolis-Hastings is available in the repo [https://github.com/compops/pmmh-correlated2015](pmmh-correlated2015). These are excellent resources for getting up to speed with the current frontier in research connected to PMH.
+
 Copyright information
 --------------
 ``` R
@@ -53,4 +57,4 @@ Copyright information
 # with this program; if not, write to the Free Software Foundation, Inc.,
 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 ##############################################################################
-```
+```

r/README.md

@@ -1,27 +1,27 @@
 # R code for PMH tutorial
 
-This R code implements the Kalman filter (KF), particle filter (PF) and particle Metropolis-Hastings (PMH) algorithm for two different dynamical models: a linear Gaussian state-space (LGSS) model and a stochastic volatilty (SV) model. Note that the Kalman filter can only be employed for the first of these two models. The details of the code is described in the tutorial paper available at: http://arxiv.org/pdf/1511.01707
+This R code implements the Kalman filter (KF), particle filter (PF) and particle Metropolis-Hastings (PMH) algorithm for two different dynamical models: a linear Gaussian state-space (LGSS) model and a stochastic volatility (SV) model. Note that the Kalman filter can only be employed for the first of these two models. The details of the code is described in the tutorial paper available at: http://arxiv.org/pdf/1511.01707
 
 Requirements
 --------------
 The code is written and tested for R 3.2.2 and makes use of the packages Quandl and mvtnorm. These can be installed in R by executing the command:
 ``` R
-install.packages(c("mvtnorm", "Quandl")) 
+install.packages(c("mvtnorm", "Quandl"))
 ```
 
 Main script files
 --------------
 These are the main script files that implement the various algorithms discussed in the tutorial.
 
-**example1-lgss.R** State estimation in a LGSS model using the KM and a fully-adapted PF (faPF). The code is discussed in Section 3.1 and the results are presented in Section 3.2 as Figure 4 and Table 1.
+**example1-lgss.R** State estimation in a LGSS model using the KF and a fully-adapted PF (faPF). The code is discussed in Section 3.1 and the results are presented in Section 3.2 as Figure 4 and Table 1.
 
 **example2-lgss.R** Parameter estimation of one parameter in the LGSS model using PMH with the faPF as the likelihood estimator. The code is discussed in Section 4.1 and the results are presented in Section 4.2 as Figure 5.
 
 **example3-sv.R** Parameter estimation of three parameters in the SV model using PMH with the bootstrap PF as the likelihood estimator. The code is discussed in Section 5.1 and the results are presented in Section 5.2 as Figure 6. The code takes about an hour to run.
 
 **example4-sv.R** Modified version of the code in *example3-sv.R* to make use of a better tailored parameter proposal. The details are discussed in Section 6.3.2 and the results are presented in the same section as Figures 7 and 8. Note that the only difference in the code is that the variable stepSize is changed.
 
-**example5-sv.R** Modified version of the code in *example3-sv.R* to make use of another parameterisation of the model and a better tailored parameter proposal. The details are discussed in Section 6.3.3 and the results are presented in the same section. Note that the differences in the code is the use of another implemenation of PMH ant that the variable stepSize is changed.
+**example5-sv.R** Modified version of the code in *example3-sv.R* to make use of another parameterisation of the model and a better tailored parameter proposal. The details are discussed in Section 6.3.3 and the results are presented in the same section. Note that the differences in the code is the use of another implementation of PMH ant that the variable stepSize is changed.
 
 
 Additional script files for creating plots for tutorial (extra-code-for-tutorial/)
@@ -55,14 +55,14 @@ These directories are placeholders for the output from running the code. The wor
 
 Adapting the code for another model
 --------------
-The code provided in *helpers/stateInference.R* and *helpers/parameterInferernce.R* is quite general. To adapt this code for your own model, you can start with the code in *example3-sv.R* together with the functions *particleFilterSVmodel* and *particleMetropolisHastings* from the helpers. 
+The code provided in *helpers/stateInference.R* and *helpers/parameterInferernce.R* is quite general. To adapt this code for your own model, you can start with the code in *example3-sv.R* together with the functions *particleFilterSVmodel* and *particleMetropolisHastings* from the helpers.
 
 ### Particle filter
 In the particle filter, you need to change the lines connected to: (i) the sampling of the initial state, (ii) the propagation of particles and (iii) the computation of the weights. For (i), you need to rewrite:
 ``` R
   particles[, 1] <- rnorm(noParticles, mu, sigmav / sqrt(1 - phi^2))
 ```
-to fit your model. Two simple choices are to make use of the stationary distribution of the state (as is done for the SV model) computed by hand or to initialise all particles to some value (as is done in the LGSS model) by:
+to fit your model. Two simple choices are to make use of the stationary distribution of the state (as is done for the SV model) computed by hand or to initialize all particles to some value (as is done in the LGSS model) by:
 ``` R
   particles[, 1] <- initialState
 ```
@@ -82,10 +82,10 @@ For (iii), you need to change:
 ```
 to something else. For the bPF, this corresponds to the observation process of your state-space model.
 
-Finally, note that the particle filter implemetation can only be used for state-space models where the state and observation are scalar. However, it is quite straightforward to make use of particle filtering when the state and/or observations are multivariate. It is basically only bookkeeping. If the dimension of the state is larger than say 5, good proposals are usually required to not run into the curse of dimensionality. This is a hot current research topic in the particle filtering literature.
+Finally, note that the particle filter implementation can only be used for state-space models where the state and observation are scalar. However, it is quite straightforward to make use of particle filtering when the state and/or observations are multivariate. It is basically only bookkeeping. If the dimension of the state is larger than say 5, good proposals are usually required to not run into the curse of dimensionality. This is a hot current research topic in the particle filtering literature.
 
 ### Particle Metropolis-Hastings
-The implemenation of the PMH algorithm is general and does not require any larger changes if the model is changed. The dimensionality of the variables *xHatFiltered*, *xHatFilteredProposed*, *theta* and *thetaProposed* needs to be altered to match the dimensionality of the state and the number of parameters in the new state-space model. Moreover, the initial value of theta and the proposal distribution need to be calibrated for your new model. The simplest way to do this is by so-called pilot runs. Set the initial value to something reasonable and stepSize to a diagonal matrix with quite small elements, so that you get at least some accepted proposed values. After the pilot run, adapt the proposal as is discussed in 6.3.2 and initialise the PMH algorithm in the estimated posterior mean. Repeat this one or two more times or until you are satisfied. 
+The implementation of the PMH algorithm is general and does not require any larger changes if the model is changed. The dimensionality of the variables *xHatFiltered*, *xHatFilteredProposed*, *theta* and *thetaProposed* needs to be altered to match the dimensionality of the state and the number of parameters in the new state-space model. Moreover, the initial value of theta and the proposal distribution need to be calibrated for your new model. The simplest way to do this is by so-called pilot runs. Set the initial value to something reasonable and stepSize to a diagonal matrix with quite small elements, so that you get at least some accepted proposed values. After the pilot run, adapt the proposal as is discussed in 6.3.2 and initialise the PMH algorithm in the estimated posterior mean. Repeat this one or two more times or until you are satisfied.
 
 It is known that this simple version of PMH performs bad when the number of parameters is larger than about 5. To circumvent this problem, see the suggestions in Sections 4.3 and 6. It is also discussed there how to choose the number of particles *noParticles* and the number of iterations *noIterations* to use in the PMH algorithm. *noBurnInIterations* can be selected by looking at the trace plot for when the Markov chain has reached its steady-state/stationarity. I usually use *noIterations* as 10,000 or 30,000 (with *noBurnInIterations* as 3,000 or 10,0000) to get good posterior estimates but these runs take time. Also, using *noParticles* as somewhere between *T* and 2*T* is a good place to start.
 
@@ -107,4 +107,4 @@ Copyright information
 # with this program; if not, write to the Free Software Foundation, Inc.,
 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 ##############################################################################
-```
+```

r/helpers/parameterEstimation.R

@@ -224,7 +224,7 @@ particleMetropolisHastingsSVmodelReparameterised <- function(y, initialTheta,
       logPrior3 <- dgamma(thetaProposed[k, 3], 3, 10, log = TRUE) - dgamma(theta[k - 1, 3], 3, 10, log = TRUE)
       logPrior <- logPrior1 + logPrior2 + logPrior3
 
-      logJacob1 <- log(abs(1 - thetaProposed[k, 2]^2)) -log(abs(1 - theta[k - 1, 2]^2))
+      logJacob1 <- log(abs(1 - thetaProposed[k, 2]^2)) - log(abs(1 - theta[k - 1, 2]^2))
       logJacob2 <- log(abs(thetaProposed[k, 3])) - log(abs(theta[k - 1, 3]))
       logJacob <- logJacob1 + logJacob2