|
nmut = int(clade.branch_length * mutation_rate * ncol) |
The current implementation uses the rounded expectation over the full sequence as the number of mutations nmut . This can be problematic in edge cases (sequences very short or mutation rate very low).
A better alternative in these cases would be drawing nmut probabilistically:
nmut = np.random.poisson(lam=clade.branch_length*mutation_rate, size=ncol).sum()