ropensci / mcbette

Model Comparison using babette
GNU General Public License v3.0
7 stars 3 forks source link

Estimate evidence for primates #25

Open richelbilderbeek opened 4 years ago

richelbilderbeek commented 4 years ago

Depends on:

richelbilderbeek commented 4 years ago

Whatever I do, I get:

epsilon particle_count marg_log_lik marg_log_lik_sd
1.000000000000000000000e+20 1 -6482.027291882671306666 1.870667062915664979883
1.000000000000000000000e+00 1 -6482.027291882671306666 1.870667062915664979883
9.999999999999999451533e-21 1 -6482.027291882671306666 1.870667062915664979883
1.000000000000000000000e+20 16 -6482.027291882671306666 1.870667062915664979883
1.000000000000000000000e+00 16 -6482.027291882671306666 1.870667062915664979883
9.999999999999999451533e-21 16 -6482.027291882671306666 1.870667062915664979883

with code:

# R script to estimate the evidence of the Primates.nex alignment
# Using different settings

library(mcbette)

options(digits = 22)

fasta_filename <- tempfile()

beastier::save_nexus_as_fasta(
  nexus_filename = beastier::get_beast2_example_filename("Primates.nex"),
  fasta_filename = fasta_filename
)

df <- expand.grid(
  epsilon = c(1e20, 1, 1e-20),
  particle_count = c(1, 16),
  marg_log_lik = NA,
  marg_log_lik_sd = NA
)

for (i in seq_len(nrow(df))) {
  print(i)
  inference_model <- create_test_ns_inference_model()
  inference_model$mcmc$epsilon <- df$epsilon[i]
  inference_model$mcmc$particle_count <- df$particle_count[i]
  inference_model$mcmc$chain_length <- 1e10
  evidence <- mcbette::est_marg_lik(
    fasta_filename = fasta_filename,
    inference_model = inference_model
  )
  df$marg_log_lik <- evidence$marg_log_lik
  df$marg_log_lik_sd <- evidence$marg_log_lik_sd
}
knitr::kable(df)
richelbilderbeek commented 4 years ago

Because every NS setting I run (for a same model, for a same alignment), I get the same marginal likelihood estimation.

Therefore, I predict this is a bug on my side.

richelbilderbeek commented 4 years ago

I created an NS Issue to ask for an example in which epsilon actually does something.

richelbilderbeek commented 4 years ago

Aha, the subchain is too short :man_facepalming:...

epsilon never has an effect :confused: