The package contains a simple Metropolis-Hastings sampler:
tab <- array(c(35L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L), dim = c(3L, 3L))
tab
#> [,1] [,2] [,3]
#> [1,] 35 0 0
#> [2,] 0 1 0
#> [3,] 0 0 0
w_b <- wgsLR::estimate_w_bayesian(tab) |> posterior_samples()
mean(w_b)
#> [1] 0.008478139Instead of using the built-in simple Metropolis-Hastings sampler, it may be a good idea to use e.g. Stan for production purposes:
library(cmdstanr)
model <- "
// https://en.wikipedia.org/wiki/Beta_distribution#Four_parameters
// https://discourse.mc-stan.org/t/four-parameter-beta-distribution/27216/2
functions {
real beta_4p_lpdf(real y, real alpha, real beta, real lwr, real upr) {
// Scale 4-parameter Beta RV to 2-parameter Beta RV
real x = (y - lwr) / (upr - lwr);
// Return scaled 2-parameter beta lpdf
return beta_lpdf(x | alpha, beta) - log(upr - lwr);
}
}
data {
int<lower = 1> N_trials;
array[3] int<lower = 0, upper = N_trials> ans;
}
parameters {
real<lower=0,upper=0.5> w;
}
transformed parameters {
simplex[3] p;
p[1] = w^4 + 4*w^2 * (1-w)^2 + (1-w)^4;
p[2] = 2*(w*(1-w))^2;
p[3] = 4 * w^3 * (1-w) + 4*w*(1-w)^3;
}
model {
w ~ beta_4p(1, 1, 0, 0.5);
ans ~ multinomial(p);
}
"
f <- write_stan_file(model)
m <- cmdstan_model(f)
ns <- wgsLR::get_ns(tab)
ns
#> [1] 36 0 0
fit_mn <- m$sample(data = list(N_trials = sum(ns), ans = ns),
chains = 4,
show_messages = FALSE)
fit_mn
#> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
#> lp__ -5.87 -5.55 0.81 0.36 -7.52 -5.29 1.00 1116 823
#> w 0.01 0.00 0.01 0.00 0.00 0.02 1.00 1094 731
#> p[1] 0.97 0.98 0.03 0.02 0.92 1.00 1.00 1094 731
#> p[2] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1094 731
#> p[3] 0.03 0.02 0.03 0.02 0.00 0.08 1.00 1094 731
fit_mn$draws("w") |> c() |> mean()
#> [1] 0.006894489Please be aware of the “Cautionary note” in the README (about not just using standard VCF-files).