### Introduction

`SAMtool`

now supports MCMC sampling of the RCM using Stan called through the `tmbstan`

R package. In this way, individual simulations within the operating model will be conditioned on individual draws of the MCMC chain.

Below is a schematic of the steps to setup the operating model using RCM and MCMC sampling. We will again use the Pacific cod example which sets up an age-structured model to mimic a delay-difference model, with knife-edge selectivity equal to knife-edge maturity.

### MPD fit

The initial fit is generated by a function call to `RCM()`

with a pre-made operating model (`OM`

) and `RCMdata`

object. This fit sets up the structure of the model, including any priors:

```
library(SAMtool)
#### Identify the maturity ogive to fix selectivity
data(pcod)
mat_ogive <- pcod$OM@cpars$Mat_age[1, , 1]
#### Remove stochastic samples of M and steepness, will estimate from prior
OM <- pcod$OM
OM@cpars$M <- OM@cpars$h <- NULL
#### Update priors for M, h, and R0 (in addition to q)
prior <- pcod$prior
prior$R0 <- c(3, exp(1), exp(12))
prior$h <- c(0.7, 0.15)
prior$M <- c(0.5, 0.1)
#### MPD fit
MPD <- RCM(OM = OM, data = pcod$data,
condition = "catch", mean_fit = TRUE,
selectivity = "free", s_selectivity = rep("SSB", ncol(pcod$data@Index)),
vul_par = matrix(mat_ogive, length(mat_ogive), 1),
map_vul_par = matrix(NA, length(mat_ogive), 1),
map_log_early_rec_dev = rep(1, OM@maxage),
prior = prior)
```

The TMB model containing the MPD fit will be in `MPD@mean_fit$obj`

.

### MCMC

Once the MPD has been optimized, the `posterior()`

function provides a convenient way to draw MCMC samples. This function is primarily a wrapper to `tmbstan::tmbstan()`

which is also a wrapper to `rstan::sampling()`

. Function arguments to `posterior()`

configure the MCMC, e.g., the number of chains, iterations, thin rate, etc., and are passed on to these two functions.

`MCMC <- posterior(MPD, chains = 2, iter = 2000, warmup = 1000, thin = 5)`

The MCMC object is a `stanfit`

object containing the parameters from the saved MCMC draws. Therefore, any methods and diagnostic functions that supports Stan models can be used with the MCMC object. These functions can be used to determine how whether the MCMC has converged.

```
#### Look at MCMC diagnostics (autocorrelation, wormplots, etc.)
sso <- shinystan::launch_shinystan(MCMC)
```

### Update RCM

Now that we have a converged MCMC, we can update the operating model from a subset of these MCMC draws with `RCMstan()`

:

```
#### Identify the simulations and chains from which to grab parameters from the MCMC object
sims <- matrix(NA, OM@nsim, 2)
set.seed(234)
sims[, 1] <- sample(MCMC@sim$chains, OM@nsim, replace = TRUE)
sims[, 2] <- sample(MCMC@sim$n_save[1], OM@nsim)
#### Update RCM output and operating model
RCM_MCMC <- RCMstan(MPD, MCMC, sims)
```

`RCMstan`

returns an updated `RCModel`

and we can generate an updated markdown report:

`plot(RCM_MCMC)`

We can now run the MSE with the conditioned operating model.

`MSE <- runMSE(RCM_MCMC@OM, ...)`