Repository: tskit-dev/msprime (language: python)

Issue:
Review and document default params for microsat models
In #2071 we added default parameters for all the new microsat models. This is convenient, but may be a bad idea as users might assume that this is *the* microsat model.

What was the reasoning for the various defaults that you chose @andrewkern?


Discussion / hints:
We also need to remove the hard-coded values from the signatures, and hopefully can migrate some of the specific values up into the superclass.
> What was the reasoning for the various defaults that you chose @andrewkern?

hopefully the final doc strings that got merged lay this out sufficiently? 

There's a few changes I'd like to make, where we're hard-coding values into the signature. I'll make a quick PR with the updates now.

Fix patch:
diff --git a/docs/mutations.md b/docs/mutations.md
index a96a06c34..bb5b4c685 100644
--- a/docs/mutations.md
+++ b/docs/mutations.md
@@ -681,19 +681,22 @@ biases.
 Consider the [Combined DNA Index System](https://en.wikipedia.org/wiki/Combined_DNA_Index_System) (CODIS), a US National database
 of genotypes maintained by the justice system for use in forensic DNA analysis.
 The CODIS system relies on 20 core microsatellite loci, that are unlinked across
-the human genome. Let's simulate a sample of 20 individuals at 20 unlinked loci,
+the human genome. Let's simulate a sample of 5 individuals at 20 unlinked loci,
 and then apply the {class}`.SMM` model to the sample.
+After running the simulation we then examine the variation across the loci, by
+iterating through the variants.
 
-Finally we will summarize variation across our loci, by
-iterating through the variants at each position / locus.
-We will print the position of the locus, then a list of allelic states, followed
-by haploid genotypes for each chromosome at that site.
+```{note}
+The alleles in {class}`.MicrosatMutationModel`s such as this represent the
+copy number of the repeat, not the actual sequences observed. See below for some basic
+techniques for working with this kind of data.
+```
 
 ```{code-cell} python
 
-ts = msprime.sim_ancestry(20, random_seed=2, sequence_length=20,
+ts = msprime.sim_ancestry(5, random_seed=2, sequence_length=20,
     recombination_rate=0.5, population_size=100_000)
-model = msprime.SMM(1, 10)
+model = msprime.SMM(lo=1, hi=10)
 mts = msprime.sim_mutations(ts, rate=1e-5, model=model, random_seed=1)
 for var in mts.variants():
     print("Site at:", var.site.position)
@@ -701,8 +704,36 @@ for var in mts.variants():
     print("    Genotypes:", var.genotypes)
 ```
 
-We can also look at the dependence of the variance in repeat
-number on mutation rate under the SMM model. Much has been
+#### Working with copy numbers
+
+In order to perform some statistical analysis of microsatellite models,
+we need a function to convert the output of msprime into something that
+we can work with easily. The variation data is returned using the
+tskit {meth}`tskit.TreeSequence.variants` method, which uses an array of
+alleles (the different states observed at a site) and an array of
+"genotypes", which are the observations at that site for all of the samples.
+Usually, we think of alleles as being nucleotides or amino acids which
+are represented as strings. In this case, however, the alleles
+represent copy-numbers which are still represented as strings. So,
+to get the observed copy numbers, we must convert the alleles
+to a numeric type, and obtain the actual observations by
+indexing into this array. This is done in the following function:
+
+```{code-cell} python
+def copy_number_matrix(ts):
+    """
+    Returns the copy number matrix from the specified tree sequence
+    simulated under a MicrosatMutationModel.
+    """
+    C = np.zeros((ts.num_sites, ts.num_samples), dtype=int)
+    for var in ts.variants():
+        alleles = np.array([int(allele) for allele in var.alleles])
+        C[var.site.id] = alleles[var.genotypes]
+    return C
+```
+
+Using this new tool, we can now look at the dependence of the variance in repeat
+number on mutation rate under the {class}`.SMM` model. Much has been
 written over the years on how to estimate the mutation
 rate from the variance in repeat number
 (e.g., [Zhivotovsky and Feldman, 1995](https://pubmed.ncbi.nlm.nih.gov/8524801/)),
@@ -711,17 +742,17 @@ so let's simulate that signal:
 ```{code-cell} python
 from matplotlib import pyplot as plt
 
-rates = [1e-2,1e-3, 1e-4, 1e-5, 1e-6]
-model = msprime.SMM()
-vars = []
+rates = np.linspace(1e-2, 1e-4, num=20)
+variances = []
 ts = msprime.sim_ancestry(100, random_seed=2, sequence_length=1, population_size=100_000)
 for r in rates:
-    mts = msprime.sim_mutations(ts, rate=r, model=model, random_seed=1)
-    vars.append(mts.genotype_matrix().var())
+    mts = msprime.sim_mutations(ts, rate=r, model=msprime.SMM(), random_seed=1)
+    C = copy_number_matrix(mts)
+    variances.append(C.var())
 
-plt.plot(rates, vars)
-plt.xlabel('mutation rate')
-plt.ylabel('variance in repeat number')
+plt.plot(rates, variances)
+plt.xlabel('Mutation rate')
+plt.ylabel('Variance in repeat number');
 ```
 
 So here we are seeing a clear relationship between the variance
diff --git a/msprime/mutations.py b/msprime/mutations.py
index 6090baab3..98e3a7fa0 100644
--- a/msprime/mutations.py
+++ b/msprime/mutations.py
@@ -20,6 +20,8 @@
 Module responsible for generating mutations on a given tree sequence.
 """
 import inspect
+import math
+import numbers
 import sys
 
 import numpy as np
@@ -220,9 +222,7 @@ def general_microsat_rate_matrix(s, u, v, p, m, lo, hi):
     :param float m: success prob of truncated gamma distribution
     :param int lo: lower bound on repeat number
     :param int hi: upper bound on repeat number
-
     :return: rate matrix
-
     """
 
     # Helper functions for the Sainudiin et al. rate matrix
@@ -307,58 +307,75 @@ class MicrosatMutationModel(MatrixMutationModel):
     """
     General class for Microsatellite mutation models
     which parameterizes the transition matrix describing
-    changes among alleles using 5 parameters according
-    to the model of Sainudiin et al. (2004).
+    changes among alleles using five parameters according
+    to the model of
+    `Sainudiin et al. (2004) <https://doi.org/10.1534/genetics.103.022665>`_
 
     Concretely, mutation rates under this model are computed as follows.
     Let ``qij`` be the transition rate from equation (3) in Sainudiin et al,
     and let ``M`` be the largest row sum of the matrix whose i,j-th entry
     is ``qij`` and whose diagonals are zero. Then if the mutation rate
     passed to :func:`.sim_mutations` is ``m``, the mutation rate from
-    allele ``i`` to allele `j`` is ``m * qij / M``.
+    allele ``i`` to allele ``j`` is ``m * qij / M``.
 
-    :param float s: strength of length dependence on mutation rate
-        (i.e. s=0 uniform rate)
-    :param float u: constant bias parameter; must be between 0 and 1
-    :param float v: linear bias parameter
-    :param float p: probability of a single step mutation
-    :param float m: 1 / (mean multistep mutation size)
-    :param int lo: lower bound on repeat number
-    :param int hi: upper bound on repeat number
+    Note that the values of ``lo`` and ``hi`` are in units of repeat number,
+    and are unrelated to repeat length
+    (e.g., the repeat itself can be dinucleotide, tri-nucleotide, etcetera).
+
+    The default values for ``lo`` and ``hi`` are chosen based on
+    FIXME.
+
+    :param float s: strength of length dependence on mutation rate.
+        Defaults to 0, i.e., no length dependence.
+    :param float u: Constant bias parameter; must be between 0 and 1.
+        Defaults to 0.5, i.e., no bias.
+    :param float v: Linear bias parameter. Defaults to 0, i.e., no bias.
+    :param float p: Probability of a single step mutation. Defaults to 1.
+    :param float m: 1 / (mean multistep mutation size). Defaults to 1.
+    :param int lo: Lower bound on repeat number. Defaults to 2.
+    :param int hi: Upper bound on repeat number (inclusive). Defaults to 50.
     :param list[float] root_distribution: An array of float values of length hi-lo+1
         which should sum to 1, and give the probabilities used to choose
         the ancestral state of each microsatellite.
         Defaults to the stationary distribution of the model.
-
     """
 
     def __init__(
         self,
-        s=0,
-        u=0.5,
-        v=0.0,
-        p=1.0,
-        m=1.0,
-        lo=2,
-        hi=50,
+        *,
+        s=None,
+        u=None,
+        v=None,
+        p=None,
+        m=None,
+        lo=None,
+        hi=None,
         root_distribution=None,
     ):
-        if np.isinf(lo):
-            raise ValueError("lo is infinite")
-        if np.isinf(hi):
-            raise ValueError("hi is infinite")
+        s = 0.0 if s is None else s
+        u = 0.5 if u is None else u
+        v = 0.0 if v is None else v
+        p = 1.0 if p is None else p
+        m = 1.0 if m is None else m
+        lo = 2 if lo is None else lo
+        hi = 50 if hi is None else hi
+
+        params = [s, u, v, p, m, lo, hi]
+        if any(not math.isfinite(x) for x in params):
+            raise ValueError("All parameters must be finite")
+        if any(not isinstance(x, numbers.Number) for x in params):
+            raise ValueError("All parameters must be numeric")
+
+        lo = int(lo)
+        hi = int(hi)
         if lo > hi:
             raise ValueError("lo > hi")
-        if not lo >= 0:
+        if lo < 0:
             raise ValueError("lo must be >= 0")
         if s <= -1 / (hi - lo + 1):
             raise ValueError("s must be > -1/(hi-lo+1)")
         if u <= 0 or u >= 1:
             raise ValueError("u must be between 0 and 1")
-        if type(lo) == float:
-            lo = int(lo)
-        if type(hi) == float:
-            hi = int(hi)
         alleles = [str(int(x)) for x in range(lo, hi + 1)]
         transition_matrix = general_microsat_rate_matrix(s, u, v, p, m, lo, hi)
         if root_distribution is None:
@@ -379,39 +396,29 @@ class SMM(MicrosatMutationModel):
     This is a :class:`.MicrosatMutationModel` with alleles ``[lo, .. , hi]``,
     a root distribution (stationary distribution by default),
     and a transition matrix that allows only mutations that
-    change the number of repeats by +/- 1: concretely, if the mutation
+    change the number of repeats by +/- 1. Concretely, if the mutation
     rate is ``m``, then an allele ``k`` mutates to ``k+1`` and ``k-1``
     at rate ``m/2`` each,
     except that mutations above ``hi`` or below ``lo`` result in no change.
 
-    Note that the values of ``lo`` and ``hi`` are in units of repeat number,
-    and are unrelated to repeat length
-    (e.g., the repeat itself can be dinucleotide, tri-nucleotide, etcetera).
-
-    :param int lo: lower bound of the microsatellite repeat copy number
-    :param int hi: upper bound of the microsatellite repeat copy number
-    :param list[float] root_distribution: An array of float values of length hi-lo+1
-        which should sum to 1, and give the probabilities used to choose
-        the ancestral state of each microsatellite.
-        Defaults to the value of stationary distribution
-        of the model.
-
+    :param int lo: Repeat number lower bound. See the documentation for
+        :class:`.MicrosatMutationModel` for details.
+    :param int hi: Repeat number upper bound. See the documentation for
+        :class:`.MicrosatMutationModel` for details.
+    :param list[float] root_distribution: See the documentation for
+        :class:`.MicrosatMutationModel` for details.
     """
 
-    def __init__(self, lo=2, hi=50, root_distribution=None):
-        # SMM specific parameters
-        s = 0
-        u = 0.5
-        v = 0
-        p = 0
-        m = 1
-        super().__init__(s, u, v, p, m, lo, hi, root_distribution)
+    def __init__(self, lo=None, hi=None, *, root_distribution=None):
+        super().__init__(
+            s=0, u=0.5, v=0, p=0, m=1, lo=lo, hi=hi, root_distribution=root_distribution
+        )
 
 
 class TPM(MicrosatMutationModel):
     """
-    Two-phase mutation model of DiRienzo et al., (1994),
-    https://doi.org/10.1073/pnas.91.8.3166
+    Two-phase mutation model of
+    `DiRienzo et al. (1994) <https://doi.org/10.1073/pnas.91.8.3166>`_
 
     This models evolution of microsatellite repeat number in a manner
     that allows for multi-step mutations in copy number of microsatellite repeats.
@@ -424,88 +431,80 @@ class TPM(MicrosatMutationModel):
     This is a :class:`.MicrosatMutationModel` with alleles ``[lo, .. , hi]``.
     For a precise definition of the mutation rates, see :class:`.MicrosatMutationModel`
     with ``s=0``, ``u=0.5``, and ``v=0``.
-    Default values of other parameters in the model
-    are those estimated in Sainudiin et al. (2004).
-
-    Note that the values of ``lo`` and ``hi`` are in units of repeat number,
-    and are unrelated to repeat length
-    (e.g., the repeat itself can be dinucleotide, tri-nucleotide, etcetera).
-
-    :param float p: probability of a single step mutation
-    :param float m: success prob of truncated gamma distribution
-    :param int lo: lower bound of the microsatellite repeat copy number
-    :param int hi: upper bound of the microsatellite repeat copy number
-    :param list[float] root_distribution: An array of float values of length hi-lo+1
-        which should sum to 1. These values are used to determine the ancestral state of
-        each microsatellite. Defaults to the stationary distribution
-        of the model
-
+    Default values for ``p`` and ``m``
+    are those estimated by
+    `Sainudiin et al. (2004) <https://doi.org/10.1534/genetics.103.022665>`_
+
+    :param float p: Probability of a single step mutation. (Default=0.9)
+    :param float m: Success prob of truncated gamma distribution. (Default=0.93)
+    :param int lo: Repeat number lower bound. See the documentation for
+        :class:`.MicrosatMutationModel` for details.
+    :param int hi: Repeat number upper bound. See the documentation for
+        :class:`.MicrosatMutationModel` for details.
+    :param list[float] root_distribution: See the documentation for
+        :class:`.MicrosatMutationModel` for details.
     """
 
-    def __init__(self, p=0.9, m=0.93, lo=2, hi=50, root_distribution=None):
+    def __init__(self, *, p=None, m=None, lo=None, hi=None, root_distribution=None):
+        p = 0.9 if p is None else p
+        m = 0.93 if m is None else m
         if not (0 < p < 1.0):
             raise ValueError("p must be between 0 and 1")
         if not (0 < m < 1.0):
             raise ValueError("m must be between 0 and 1")
-        # TPM-specific fixed params
-        s = 0
-        u = 0.5
-        v = 0
-        super().__init__(s, u, v, p, m, lo, hi, root_distribution)
+        super().__init__(
+            s=0, u=0.5, v=0, p=p, m=m, lo=lo, hi=hi, root_distribution=root_distribution
+        )
 
 
 class EL2(MicrosatMutationModel):
     """
-    Equal rate, linear biased, two-phase mutation model of Garza et al., (1995),
-    [10.1093/oxfordjournals.molbev.a040239] as parameterized by
-    the Sainudiin et al. (2004) model.
+    Equal rate, linear biased, two-phase mutation model of
+    `Garza et al. (1995) <10.1093/oxfordjournals.molbev.a040239>`_
+    as parameterized by
+    `Sainudiin et al. (2004) <https://doi.org/10.1534/genetics.103.022665>`_
 
     This models evolution of microsatellite repeat number in a manner
-    that allows for multi-step mutations in copy number of microsatellite repeats.
+    that allows for multi-step mutations in copy number.
     This is parameterized by `p` the probability of a single step mutation,
     and `m` the success probability of the truncated gamma distribution
-    describing the distribution of longer steps such that the
+    describing the distribution of longer steps, such that the
     mean multi-step mutation is length `1/m`.
     For this model both `p` and `m` need to be set to < 1.0.
 
-    In addition EL2 has a constant bias parameter `u`, where
-    :math:`0 < u < 1` and a linear bias parameter `v`, where
-    :math:`-\\infty < v < \\infty`.
+    See :class:`.MicrosatMutationModel` for details on the  constant
+    bias parameter `u` and a linear bias parameter `v`.
 
     This is a :class:`.MicrosatMutationModel` with alleles ``[lo, .. , hi]``
-    For a precise definition of the mutation rates, see :class:`.MicrosatMutationModel`
-    with ``s=0``.
-    Default values for other parameters of the model
-    are those estimated in Sainudiin et al. (2004).
-
-    Note that the values of ``lo`` and ``hi`` are in units of repeat number,
-    and are unrelated to repeat length
-    (e.g., the repeat itself can be dinucleotide, tri-nucleotide, etcetera).
-
-    :param float p: probability of a single step mutation
-    :param float m: success prob of truncated gamma distribution
-    :param float u: constant bias parameter
-    :param float v: linear bias parameter
-    :param int lo: lower bound of the microsatellite repeat copy number
-    :param int hi: upper bound of the microsatellite repeat copy number
-    :param list[float] root_distribution: An array of float values of length hi-lo+1
-        which should sum to 1. These values are used to determine the ancestral state of
-        each microsatellite. Defaults to the stationary distribution
-        of the model
-
+    For a precise definition of the mutation rates,
+    see :class:`.MicrosatMutationModel`
+    with ``s=0``. Default values for ``m``, ``u`` and ``v``
+    are those estimated by Sainudiin et al. (2004).
+
+    :param float m: Success prob of truncated gamma distribution (0.43).
+    :param float u: Constant bias parameter (Default=0.68).
+    :param float v: Linear bias parameter (Default=0.037).
+    :param int lo: Repeat number lower bound. See the documentation for
+        :class:`.MicrosatMutationModel` for details.
+    :param int hi: Repeat number upper bound. See the documentation for
+        :class:`.MicrosatMutationModel` for details.
+    :param list[float] root_distribution: See the documentation for
+        :class:`.MicrosatMutationModel` for details.
     """
 
-    def __init__(self, m=0.43, u=0.68, v=0.037, lo=2, hi=50, root_distribution=None):
+    def __init__(
+        self, *, m=None, u=None, v=None, lo=None, hi=None, root_distribution=None
+    ):
+        m = 0.43 if m is None else m
+        u = 0.68 if u is None else u
+        v = 0.037 if v is None else v
         if not (0 < m < 1.0):
             raise ValueError("m must be between 0 and 1")
         if not (0 < u < 1.0):
             raise ValueError("u must be between 0 and 1")
-        if not (-np.inf < v < np.inf):
-            raise ValueError("v must be between -inf and inf")
-        # EL2 specific params
-        p = 0.0
-        s = 0
-        super().__init__(s, u, v, p, m, lo, hi, root_distribution)
+        super().__init__(
+            s=0, u=u, v=v, p=0, m=m, lo=lo, hi=hi, root_distribution=root_distribution
+        )
 
 
 class JC69(MatrixMutationModel):


Test patch:
diff --git a/tests/test_mutations.py b/tests/test_mutations.py
index ab8ddfb47..36326ccd8 100644
--- a/tests/test_mutations.py
+++ b/tests/test_mutations.py
@@ -412,6 +412,22 @@ def test_EL2(self):
                 exp = sum_ij2 / sum_ij1
                 assert np.abs(exp - a_ratio) < 1e-6
 
+    @pytest.mark.parametrize(
+        "m, u, v, p, s, lo, hi",
+        [
+            (1.0, 0.62, 0.015, 0, 0.88, 1, 10.0),
+            (0.43, 0.68, 0.037, 1.0, 0.5, 3.0, 10),
+            (0.83, 0.58, 0.37, 0, 0.88, 1, 20),
+            (0.03, 0.18, 0.5, 1.0, 1, 1, 4),
+        ],
+    )
+    def test_MicrosatMutationModel(self, m, u, v, p, s, lo, hi):
+        model = msprime.MicrosatMutationModel(s=s, u=u, v=v, p=p, m=m, lo=lo, hi=hi)
+        self.validate_model(model)
+        self.validate_stationary(model)
+
+
+class TestMicrosatModels:
     @pytest.mark.parametrize(
         "m, u, v, p, s",
         [
@@ -421,10 +437,10 @@ def test_EL2(self):
             (0.43, 0.68, 0.037, 1.0, 0.88),
         ],
     )
-    def test_MicrosatMutationModel(self, m, u, v, p, s):
+    def test_transition_matrix(self, m, u, v, p, s):
         lo = 1
         hi = 10
-        model = msprime.MicrosatMutationModel(s, u, v, p, m, lo, hi)
+        model = msprime.MicrosatMutationModel(s=s, u=u, v=v, p=p, m=m, lo=lo, hi=hi)
         row_sums = model.transition_matrix.sum(axis=1)
         d = row_sums - model.transition_matrix.diagonal()
         betas = [1 + ((i - lo) * s) for i in range(lo, hi + 1)]
@@ -447,41 +463,30 @@ def test_MicrosatMutationModel(self, m, u, v, p, s):
             (0.43, 10, 0.037, 0, 0.88, 1, 5),
         ],
     )
-    def test_bad_MicrosatMutationModel(self, m, u, v, p, s, lo, hi):
+    def test_bad_params(self, m, u, v, p, s, lo, hi):
         with pytest.raises(ValueError):
-            msprime.MicrosatMutationModel(s, u, v, p, m, lo, hi)
+            msprime.MicrosatMutationModel(s=s, u=u, v=v, p=p, m=m, lo=lo, hi=hi)
 
-    @pytest.mark.parametrize(
-        "m, u, v, p, s, lo, hi",
-        [
-            (1.0, 0.62, 0.015, 0, 0.88, 1, 10.0),
-            (0.43, 0.68, 0.037, 1.0, 0.5, 3.0, 10),
-            (0.83, 0.58, 0.37, 0, 0.88, 1, 20),
-            (0.03, 0.18, 0.5, 1.0, 1, 1, 4),
-        ],
-    )
-    def test_good_params_MicrosatMutationModel(self, m, u, v, p, s, lo, hi):
-        model = msprime.MicrosatMutationModel(s, u, v, p, m, lo, hi)
-        self.validate_model(model)
-        self.validate_stationary(model)
+    def test_bad_TPM_p_value(self):
+        with pytest.raises(ValueError, match="p must be"):
+            msprime.TPM(p=1.1)
 
-    def test_bad_TPM_params(self):
-        with pytest.raises(ValueError):
-            msprime.TPM(1.1, 0.9, 2, 9)
-        with pytest.raises(ValueError):
-            msprime.TPM(0.9, 1.1, 2, 9)
+    def test_bad_TPM_m_value(self):
+        with pytest.raises(ValueError, match="m must be"):
+            msprime.TPM(m=1.1)
 
-    def test_bad_EL2_params(self):
-        with pytest.raises(ValueError):
-            msprime.EL2(1.1, 0.9, 0.1, 0.1, 2, 9)
-        with pytest.raises(ValueError):
-            msprime.EL2(0.9, 1.1, 0.1, 0.1, 2, 9)
-        with pytest.raises(ValueError):
-            msprime.EL2(0.9, 0.9, -0.1, 0.1, 2, 9)
-        with pytest.raises(ValueError):
-            msprime.EL2(0.9, 0.9, 0.1, np.inf, 2, 9)
-        with pytest.raises(ValueError):
-            msprime.EL2(0.9, 0.9, np.inf, np.inf, 2, 9)
+    def test_bad_EL2_m_value(self):
+        with pytest.raises(ValueError, match="m must be"):
+            msprime.EL2(m=1.1)
+
+    def test_bad_EL2_u_value(self):
+        with pytest.raises(ValueError, match="u must be"):
+            msprime.EL2(u=1.1)
+
+    @pytest.mark.parametrize("v", [np.inf, np.nan])
+    def test_bad_EL2_v_value(self, v):
+        with pytest.raises(ValueError, match="finite"):
+            msprime.EL2(v=v)
 
 
 class TestMutateRateMap: