Pairwise issue 388

[osorensen]

Here is "proof" that it reproduces Metropolis-Hastings output. Note that we need quite a lot of MCMC steps.

library(BayesMallows)
library(patchwork)
set.seed(1)

dat <- subset(beach_preferences, assessor < 5)

mod_init <- compute_mallows(
  data = setup_rank_data(preferences = dat),
  compute_options = set_compute_options(nmc = 10000, burnin = 1000)
)
alpha <- numeric()

mod <- mod_init
for(i in 6:60) {
  print(i)
  mod <- update_mallows(
    model = mod,
    new_data = setup_rank_data(preferences = subset(beach_preferences, assessor == i),
                               timepoint = i),
    smc_options = set_smc_options(
      n_particles = 10000, mcmc_steps = 200)
  )
  alpha <- c(alpha, mean(mod$alpha_samples))
}
#> [1] 6
#> [1] 7
#> [1] 8
#> [1] 9
#> [1] 10
#> [1] 11
#> [1] 12
#> [1] 13
#> [1] 14
#> [1] 15
#> [1] 16
#> [1] 17
#> [1] 18
#> [1] 19
#> [1] 20
#> [1] 21
#> [1] 22
#> [1] 23
#> [1] 24
#> [1] 25
#> [1] 26
#> [1] 27
#> [1] 28
#> [1] 29
#> [1] 30
#> [1] 31
#> [1] 32
#> [1] 33
#> [1] 34
#> [1] 35
#> [1] 36
#> [1] 37
#> [1] 38
#> [1] 39
#> [1] 40
#> [1] 41
#> [1] 42
#> [1] 43
#> [1] 44
#> [1] 45
#> [1] 46
#> [1] 47
#> [1] 48
#> [1] 49
#> [1] 50
#> [1] 51
#> [1] 52
#> [1] 53
#> [1] 54
#> [1] 55
#> [1] 56
#> [1] 57
#> [1] 58
#> [1] 59
#> [1] 60

mod_bmm <- compute_mallows(
  data = setup_rank_data(preferences = beach_preferences),
  compute_options = set_compute_options(nmc = 50000, burnin = 1000)
)

plot(mod_bmm) + plot(mod) + plot_layout(ncol = 1)

^{Created on 2024-03-06 with reprex v2.1.0}