Skip to content

Commit

Permalink
Update comments in get_transition_probs
Browse files Browse the repository at this point in the history
  • Loading branch information
szhan committed Feb 27, 2024
1 parent 1371779 commit e5f5e3d
Showing 1 changed file with 6 additions and 6 deletions.
12 changes: 6 additions & 6 deletions python/tests/beagle_numba.py
Original file line number Diff line number Diff line change
Expand Up @@ -280,32 +280,32 @@ def get_mismatch_probs(pos, error_rate):

def get_transition_probs(cm, h, ne):
"""
Compute transition probabilities at genotyped positions.
Compute probabilities of transitioning to a random state at genotyped positions.
Note that in BEAGLE 4.1, the default value of `ne` is set to 1e6,
whereas in BEAGLE 5.4, the default value of `ne` is set to 1e5.
In BEAGLE, this default value was optimized empirically on
In BEAGLE 4.1/5.4, this default value was optimized empirically on
datasets from large outbred human populations. Also, in IMPUTE5,
the default value of `ne` is set to 1e6.
If `h` is small and `ne` is large, the transition probabilities are ~1.
In such cases, it may be desirable to set `ne` to a small value
to avoid switching between haplotypes all the time.
When using `_tskit.LsHmm`, this is `rho`.
This corresponds to `rho` in `_tskit.LsHmm`.
:param numpy.ndarray cm: Genetic map positions (cM).
:param int h: Number of reference haplotypes.
:param float ne: Effective population size.
:return: Transition probabilities.
:rtype: numpy.ndarray
"""
# Recombination rate of first site is always 0.
# E(number of crossover events) at first site is always 0.
r = np.zeros(len(cm), dtype=np.float64)
r[1:] = np.diff(cm)
c = -0.04 * (ne / h)
rho = -1 * np.expm1(c * r)
return rho
trans_probs = -1 * np.expm1(c * r)
return trans_probs


@njit
Expand Down

0 comments on commit e5f5e3d

Please sign in to comment.