From 13afd05de9052a39f271c418a377c15d529a8b54 Mon Sep 17 00:00:00 2001 From: szhan Date: Tue, 18 Jun 2024 11:53:30 +0100 Subject: [PATCH] Check that the number of haplotypes meets a minimum --- lshmm/core.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/lshmm/core.py b/lshmm/core.py index 47f6c16..c3d65b8 100644 --- a/lshmm/core.py +++ b/lshmm/core.py @@ -409,6 +409,9 @@ def get_index_in_emission_matrix_diploid(ref_genotype, query_genotype): # Miscellaneous functions. def estimate_mutation_probability(num_haps): """Return the mutation probability as defined by A2 and A3 in Li & Stephens (2003).""" + if num_haps < 3: + err_msg = "Number of haplotypes must be at least 3." + raise ValueError(err_msg) theta_tilde = 1 / np.sum([1 / k for k in range(1, num_haps - 1)]) prob_mutation = 0.5 * (theta_tilde / (num_haps + theta_tilde)) return prob_mutation