From b3c9f6b367ead408a64c4d1cd27a4781c80292f4 Mon Sep 17 00:00:00 2001 From: "Aaron A. King" Date: Mon, 7 Aug 2023 15:49:26 -0400 Subject: [PATCH] pomp v5.4; hindcasting in getting started --- _posts/2023-08-07-pomp-version-5.4.md | 17 ++++++++++++++ vignettes/getting_started.Rmd | 34 +++++++++++++++++++++++++++ vignettes/getting_started.html | 33 ++++++++++++++++---------- 3 files changed, 72 insertions(+), 12 deletions(-) create mode 100644 _posts/2023-08-07-pomp-version-5.4.md diff --git a/_posts/2023-08-07-pomp-version-5.4.md b/_posts/2023-08-07-pomp-version-5.4.md new file mode 100644 index 00000000..2e6b6b0a --- /dev/null +++ b/_posts/2023-08-07-pomp-version-5.4.md @@ -0,0 +1,17 @@ +--- +date: 7 August 2023 +layout: pomp +title: pomp version 5.4 released +--- + +**pomp** version 5.4 has been released to CRAN and is making its way to [a mirror near you](https://cran.r-project.org/mirrors.html). + +This release contains just one feature enhancement, but it eliminates an annoying warning message on some platforms that inadvertently slipped into the last round of revisions. + +### Feature enhancements + +The archiving function `stew` now saves timing information by default. + +### Bug fixes + +A bug resulting in a spurious warning about conversion of `NULL` elements that appears on some platforms has been fixed. diff --git a/vignettes/getting_started.Rmd b/vignettes/getting_started.Rmd index 58d53b81..f2613a35 100644 --- a/vignettes/getting_started.Rmd +++ b/vignettes/getting_started.Rmd @@ -1556,6 +1556,23 @@ Comment on your findings. ## Hindcast and smoothing +To this point, we have focused on the particle filter as a means of computing the likelihood, for use in estimation of unknown parameters. +It also gives us a way of estimating the likely trajectories of the latent state process. +The `filter_mean`, `pred_mean`, and `pred_var` functions give access to the *filtering* and *prediction* distributions of the latent state at each observation time. +The prediction distribution is that of $X_n\;\vert\;Y_1,\dots,Y_{n-1}$, i.e., the distribution of the latent state *conditional on* the observations at all times up to *but not including* $t_n$. +The filter distribution, by contrast, incorporates the information of the $n$-th observation: +it is the distribution of $X_n\;\vert\;Y_1,\dots,Y_n$. + +One can think of the filtering distribution as a kind of *hindcast*, i.e., an estimate of the state of the latent process at earlier times. +Although it is computed essentially at no additional cost, it is typically of limited value as a hindcast, since the amount of information used in the estimate decreases as one goes farther back into the past. +A better hindcast is afforded by the *smoothing distribution*, which is the distribution of $X_n\;\vert\;Y_1,\dots,Y_N$, i.e., that of the latent state conditional on *all* the data. + +Computation of the smoothing distribution requires more work. +Specifically, a single particle chosen at random from a particle fitler, with its entire history, is a draw from the sampling distribution [@Andrieu2010]. +To build up a picture of the sampling distribution, therefore, one can run some number of independent particle filters, sampling a single particle's trajectory from each. +The following codes accomplish this. +Note the use of the `filter_traj` function, which extracts the trajectory of a single, randomly chosen, particle. + ```{r hindcast1} r_prof |> group_by(r) |> @@ -1611,6 +1628,14 @@ quants1 |> labs(y="N") ``` +Note also that the individual trajectories must we weighted by their Monte Carlo likelihoods (see the use of `wquant` in the above). + +The plot above shows the uncertainty surrounding the hindcast. +Importantly, this includes uncertainty due to both measurement error and process noise. +It does not, however, incorporate parametric uncertainty. +To fold this additional uncertainty into the estimates requires further work. +We discuss this [in a later section](#hindcast-and-smoothing-with-parametric-uncertainty). + ## Sampling the posterior using particle Markov chain Monte Carlo ### Running multiple `pmcmc` chains @@ -1768,6 +1793,15 @@ traces |> ### Hindcast and smoothing with parametric uncertainty +Let us revisit the matter of hindcasting. +In a [preceding section](#hindcast-and-smoothing), using an ensemble of particle filters, we were able to estimate the likely trajectory that the latent state process followed in generating the data we see, *conditional on a given parameter vector*. +Naturally, since our estimate of the parameters is itself uncertain, our actual uncertainty regarding the hindcast is somewhat larger. +We can use particle Markov chain Monte Carlo to fold in the parametric uncertainty. +The following codes select one randomly-chosen particle trajectory from each iteration of our pMCMC chains. +The chains are thinned and the burn-in period is discarded, of course. +Note that it is not necessary to weight the samples by the likelihood: +if the pMCMC chains have converged, the samples can be taken to be equally weighted samples from the smoothing distribution. + ```{r hindcast2} chains |> filter_traj() |> diff --git a/vignettes/getting_started.html b/vignettes/getting_started.html index 1e8ddc7b..c8b89793 100644 --- a/vignettes/getting_started.html +++ b/vignettes/getting_started.html @@ -1709,7 +1709,7 @@

Structure of a POMP

These conditions can be represented schematically by the following diagram, which indicates the direct dependencies among model variables.

-**Conditional independence graph of a POMP.** The latent state $X_n$ at time $t_n$ is conditionally independent of its history given $X_{n-1}$. The observation $Y_n$ is conditionally independent of all other variables given $X_n$. The underlying time scale can be taken to be either discrete or continuous and the observation times need not be equally spaced. +**Conditional independence graph of a POMP.** The latent state $X_n$ at time $t_n$ is conditionally independent of its history given $X_{n-1}$. The observation $Y_n$ is conditionally independent of all other variables given $X_n$. The underlying time scale can be taken to be either discrete or continuous and the observation times need not be equally spaced.

Conditional independence graph of a POMP. The latent state \(X_n\) at time \(t_n\) is conditionally independent of its history given \(X_{n-1}\). The observation \(Y_n\) is conditionally independent of all other variables given \(X_n\). The underlying time scale can be taken to be either discrete or continuous and the observation times need not be equally spaced.

@@ -2765,7 +2765,7 @@

Estimating the log likelihood

Now, although we have observed the intended improvement in the log likelihood, we should be careful to note that the log likelihood displayed in this plot is the log likelihood of the perturbed model. This model differs from the one we are interested in. To compute the likelihood of our focal model at the parameter returned by mif2, we need to perform a few particle filter operations:

vpM |> pfilter() |> logLik() |> replicate(n=5) |> logmeanexp(se=TRUE,ess=TRUE)
##          est           se          ess 
-## -143.8135764    0.1290785    4.6997258
+## -143.9880041 0.1035849 4.7923226

Modern computers have multiple processors (cores). To take advantage of these, we can parallelize the above particle-filter computation using the circumstance and doFuture packages:

library(circumstance)
 library(doFuture)
@@ -2773,7 +2773,7 @@ 
Estimating the log likelihood

 
 vpM |> pfilter(Nrep=5) |> logLik() |> logmeanexp(se=TRUE,ess=TRUE)
##          est           se          ess 
-## -143.8173736    0.1185375    4.7326512
+## -144.0370140 0.1467573 4.6030137

In the above, the pfilter function is provided by circumstance. It causes the particle filters to be run in parallel, each using a different core. Internally, circumstance uses the doFuture package to parallelize computations. doFuture provides a number of other ways of parallelizing computations, inclusing multisession (needed on Windows computers) and cluster (for parallelizing across multiple machines). These are set using the plan() function. doFuture also provides parallel random-number generators, so that we can consider each of the parallel computations to be independent.

@@ -2977,7 +2977,7 @@

Likelihood profile

starts |>
   select(r,sigma,K,N_0,b) |>
   plot_matrix()
-

+

vpM |>
   mif2(
     starts=starts,
@@ -3031,7 +3031,7 @@ 
Likelihood profile

   geom_smooth(method="loess",span=0.2)+
   scale_x_log10()+
   labs(y=expression(sigma))
-

+

Exercise

@@ -3042,6 +3042,9 @@

Exercise

Hindcast and smoothing

+

To this point, we have focused on the particle filter as a means of computing the likelihood, for use in estimation of unknown parameters. It also gives us a way of estimating the likely trajectories of the latent state process. The filter_mean, pred_mean, and pred_var functions give access to the filtering and prediction distributions of the latent state at each observation time. The prediction distribution is that of \(X_n\;\vert\;Y_1,\dots,Y_{n-1}\), i.e., the distribution of the latent state conditional on the observations at all times up to but not including \(t_n\). The filter distribution, by contrast, incorporates the information of the \(n\)-th observation: it is the distribution of \(X_n\;\vert\;Y_1,\dots,Y_n\).

+

One can think of the filtering distribution as a kind of hindcast, i.e., an estimate of the state of the latent process at earlier times. Although it is computed essentially at no additional cost, it is typically of limited value as a hindcast, since the amount of information used in the estimate decreases as one goes farther back into the past. A better hindcast is afforded by the smoothing distribution, which is the distribution of \(X_n\;\vert\;Y_1,\dots,Y_N\), i.e., that of the latent state conditional on all the data.

+

Computation of the smoothing distribution requires more work. Specifically, a single particle chosen at random from a particle fitler, with its entire history, is a draw from the sampling distribution (Andrieu, Doucet, and Holenstein 2010). To build up a picture of the sampling distribution, therefore, one can run some number of independent particle filters, sampling a single particle’s trajectory from each. The following codes accomplish this. Note the use of the filter_traj function, which extracts the trajectory of a single, randomly chosen, particle.

r_prof |>
   group_by(r) |>
   filter(rank(-loglik)<=2) |>
@@ -3096,6 +3099,8 @@ 
Hindcast and smoothing

   geom_point(data=parus,aes(x=year,y=pop))+
   labs(y="N")

+

Note also that the individual trajectories must we weighted by their Monte Carlo likelihoods (see the use of wquant in the above).

+

The plot above shows the uncertainty surrounding the hindcast. Importantly, this includes uncertainty due to both measurement error and process noise. It does not, however, incorporate parametric uncertainty. To fold this additional uncertainty into the estimates requires further work. We discuss this in a later section.

Sampling the posterior using particle Markov chain Monte Carlo

@@ -3322,6 +3327,7 @@

Examining the posterior density

Hindcast and smoothing with parametric uncertainty

+

Let us revisit the matter of hindcasting. In a preceding section, using an ensemble of particle filters, we were able to estimate the likely trajectory that the latent state process followed in generating the data we see, conditional on a given parameter vector. Naturally, since our estimate of the parameters is itself uncertain, our actual uncertainty regarding the hindcast is somewhat larger. We can use particle Markov chain Monte Carlo to fold in the parametric uncertainty. The following codes select one randomly-chosen particle trajectory from each iteration of our pMCMC chains. The chains are thinned and the burn-in period is discarded, of course. Note that it is not necessary to weight the samples by the likelihood: if the pMCMC chains have converged, the samples can be taken to be equally weighted samples from the smoothing distribution.

chains |>
   filter_traj() |>
   melt() |>
@@ -3372,7 +3378,7 @@ 
Simulation of the fitted model

   labs(alpha="")+
   geom_line()+
   theme_bw()
-

+

The first lines above simply extract the maximum likelihood estimates (mle) from our profle computation. The next pair of lines plug these MLE parameters into a ‘pomp’ object (mlepomp) containing the model and the data. The last set of lines do the simulation and the plotting.

Although it is clear from these plots that the estimated model has more variability and is thus able to explain the data better, it can be hard to read much from spaghetti plots such as this. It’s almost always a good idea to plot the data together with several simulated realizations in order to help assess how similar the two are.

mlepomp |>
@@ -3384,7 +3390,7 @@ 
Simulation of the fitted model

   geom_line()+
   facet_wrap(~.id)+
   theme_bw()
-

+

Model checking with probes

@@ -3411,17 +3417,17 @@

Model checking with probes

## ## $quantiles ## mean q.5% q.25% q.50% q.75% q.95% acf[1] acf[3] -## 0.495 0.415 0.605 0.280 0.375 0.825 0.695 0.825 +## 0.430 0.425 0.630 0.220 0.265 0.790 0.665 0.815 ## ## $pvals ## mean q.5% q.25% q.50% q.75% q.95% acf[1] acf[3] -## 0.9950249 0.8358209 0.7562189 0.5671642 0.7562189 0.3582090 0.6169154 0.3582090 +## 0.8656716 0.8557214 0.7263682 0.4477612 0.5373134 0.4278607 0.6766169 0.3781095 ## ## $synth.loglik -## [1] -18.70788 +## [1] -19.11099

Evidently, summary returns a list with several elements. The quantiles element contains, for each probe, what fraction of the nsim simulations had probe values below the value of the probe applied to the data. The pvals element contains \(P\)-values associated with the two-sided test of the hypothesis that the data were generated by the model.

plot(vp_probe)
-

+

The plot depicts the multivariate distribution of the probes under the model, with the data-values superimposed. On the diagonal, we see the marginal distributions of the individual probes, represented as histograms, with the vertical line signifying the value of the corresponding probe on the data. Above the diagonal, the scatterplots show the pairwise distributions of probes and the crosshairs, the corresponding data-values. Below the diagonal, the panels contains the pairwise correlations among the simulated probes.

@@ -3443,7 +3449,7 @@

Session information

pomp -5.3.1.0 +5.3.1.3 circumstance @@ -3577,6 +3583,9 @@

Session information

References

+
+

Andrieu C, Doucet A, Holenstein R (2010). “Particle Markov Chain Monte Carlo Methods.” J R Stat Soc B, 72(3), 269–342. https://doi.org/10.1111/j.1467-9868.2009.00736.x.

+

Arulampalam MS, Maskell S, Gordon N, Clapp T (2002). “A Tutorial on Particle Filters for Online Nonlinear, Non-Gaussian Bayesian Tracking.” IEEE Trans Signal Process, 50, 174–188. https://doi.org/10.1109/78.978374.