nan uncertainty estimates

[proteneer]

For certain types of overlapping distributions, the uncertainty estimate returned by BAR can benan

%matplotlib inline
import numpy as np
import pymbar
from matplotlib import pyplot as plt
for _ in range(10):
    fwd = np.random.normal(0, 10.0, size=50000)
    rev = np.random.normal(0, 175.0, size=50000)
    print(pymbar.BAR(fwd, rev))

Returns

(0.03461197183239051, nan)
(0.01789748600329588, nan)
(-0.041729050167845116, 0.006057615238226341)
(-0.0008426040270812507, nan)
/home/yutong/venv37/lib/python3.7/site-packages/numpy/core/_methods.py:48: RuntimeWarning: overflow encountered in reduce
  return umr_sum(a, axis, dtype, out, keepdims, initial, where)
(0.011851492025385824, nan)
(-0.011125827676117694, nan)
(-0.0011024556901020333, 0.006058115351639581)
(-0.14948909238990282, 0.006102804127491066)
(0.03077049043731961, 0.006088478465760845)
(-0.10509980343390524, nan)

The distributions look like:

Edit: G0 and G1 below is wrong, confused the distributions, there is probably no issue with the bias!

Futhermore, since the analytical partition functions are Z(sig)=1/(sig*sqrt(2pi)), so using G=-ln(Z) the analytical free energies are G0=3.221 and G1=6.083, The estimator also appears to be biased.

[maxentile]

Quick question: Can you confirm if this error can also be triggered using the gaussian_work test system?

pymbar/pymbar/testsystems/gaussian_work.py

Line 4 in 15f932a

def gaussian_work_example(N_F=200, N_R=200, mu_F=2.0, DeltaF=None, sigma_F=1.0, seed=None):

This would be helpful in investigating the bias part of this issue, since I'm currently unclear on how G0 and G1 are derived in the case where forward and reverse work mean and stddev all vary independently.

---

(If the error cannot be reproduced using the gaussian_work test system, it may be useful to expand this issue with a reproducing example in a different format. (E.g. by providing a pair of reduced potential energy functions u_A and u_B, and accompanying samples from p_A and p_B.))

[proteneer]

Oops yeah the estimates for G0 and G1 are clearly wrong here. Let me try to pare down from the actual distributions and not the work values.

[jchodera]

I've just suggested we also add a method to compute the number of effective samples: #427

I'm guessing these cases are seeing a collapse in the number of effective samples, meaning the estimates will be unreliable.

We should add some protection to make sure we don't return nans, thought. And we may want to warn when there is a sample size collapse.

[mrshirts]

@jchodera - see the discussion above - the example was using U_1 and U_2, not U_1-U2 and U_2-U1, so the results were nonsensical because the data was not physical. Good question about what should be returned in those cases, though, and if they should be trapped (if possible? Might be hard to distinguish in many cases from numerically bad data).

[jchodera]

Oh! Thanks for pointing that out, @mrshirts.

Could we do a quick test of consistency with the CFT? If there is sufficient data (like there was here), it can quickly flag that there's an input issue.

[mrshirts]

Could we do a quick test of consistency with the CFT? If there is sufficient data (like there was here), it can quickly flag that there's an input issue.

That's a great question. Should be able to take the log of the distribution ratio and get a straight line, and then can test how far off it is. But the tolerance is a real question - HOW far off can it be? That will take some experimentation which will take a little time.

[jchodera]

Should be able to take the log of the distribution ratio and get a straight line, and then can test how far off it is.

You're referring to tests like this paper, which plots log[p_r(-w) / p_f(+w)] and attempts to test whether this is a straight line consistent with the estimated free energy difference:

This could work, but only with a ton of samples, since histogram or kernel estimators for p(w) are going to be very noisy.

Lower-order moments of p(w) are likely to be easier to test. We could instead use the relationship

p_f(w) = p_r(-w) exp[w - Δf]

which can be expressed in terms of arbitrary functions of the work g(w):

\int dw f(w) p_f(w) = \int dw f(w) p_r(-w) exp[w - Δf]
E_f[g(w)] = E_r[g(-w) exp[Δf - w]]

to give us a relationship between moments n = 1, 2, ...:

E_f[w^n] = E_r[(-w)^n exp[Δf - w]]

(We may want to use central moments, (w - Δf)^n, or other similar functions instead.)

These should be much better conditioned and easier to check for statistical violation.

We could truncate this sum at a number of moments that are sensible based on the quantity of data and signal a warning if things are clearly in violation.

That should catch issues like this one.

• Note that I may have messed up some signs.

[mrshirts]

This could work, but only with a ton of samples, since histogram or kernel estimators for p(w) are going to be very noisy.

You can do maximum likelihood estimates of the line parameters, you don't need kernel density or histogram - see the equivalent derivation for physical_validation: https://pubs.acs.org/doi/abs/10.1021/ct300688p.

[jchodera]

I'm not entirely sure how that works with work values (where the forms of p_f(w) and p_r(w) are unknown) and arbitrary generalized reduced free energy differences. Is there a trick to deriving something like this log likelihood function only in terms of work values for this case?

[mrshirts]

The probabilities are unknown, but the ratio is known (rather, known slope but unknown constant), so it works out. In ensemble validation

ln P(U_1)/P(U_2) = ln [Q_1^-1 (\beta) Omega_1 exp(-beta_1 U) / Q_1^-2(beta_1) Omega_2 exp(-beta_2 U) 
                            = ln (Q_2/Q_1) - (beta _1 - beta_2) U 

Crooks is:

ln P(W_f)/P(-W_r) = -f + beta W_f 

So same math should work. I can get this done for 4.0 release (sometime this semester? A couple of big deadlines first).