harrispopgen / mushi

[mu]tation [s]pectrum [h]istory [i]nference
https://harrispopgen.github.io/mushi/
MIT License
24 stars 5 forks source link

Eta fit #32

Closed wsdewitt closed 4 years ago

wsdewitt commented 4 years ago

This PR implements fitting of eta(t) assuming a constant mutation rate using the method mushi.kSFS.infer_η(), and eliminates dependence on dadi. The idea is to infer eta(t) first, then condition on this eta(t) for our inference of mu(t).

Inference of eta(t) for 1KG CEU gives results roughly consistent with the literature: image

We achieve a much better fit to the total SFS than dadi, so the inferred mutation spectrum history has less artifactual scale problems (but still we see some, due the assumption of constant total mu for the eta(t) fitting step): image

The implementation relies on spicy.minimize with numerical differentiation for the eta(t) fitting step, so it's rather slow. We'll want to implement analytical gradients to speed this up.