I'm trying to test whether a likelihood has converged to walk around a single number or if it is increasing. In my (EM) algorithm, the noiseless likelihood cannot decrease, only ever increase. In my case, the likelihood is estimated and so there is noise. Therefore, each likelihood has an error on it, and therefore it can decrease at points due to noise. Also, it would be nice to detect significant decreases, as this means there is a bug in the algorithm.

So I need to test to see if there is a gradient in chain:

$P(-t \leq m \leq t) > p$,

where $t$ is the tolerance (what changes in the slope $m$ do I care about) and $p$ is the probability of convergence.

I thought the best way to do this was to fit a straight line and perform a t-test with a threshold.

Below is a MWE for testing for a non-existent gradient with data generated from a normal distribution with each data point having an uncertainty of 0.01. I also set the tolerance $t = 0.01$.

The pvalue is then estimated as 0, so I'm doing something wrong here. My question is: Is my mistake a theoretical or coding one? Is this the right way to do this?



from scipy.stats import t as studentT
from scipy.optimize import curve_fit
import numpy as np

def studentt_test(p, perr, tolerance, n):
    Returns p value for a measured p +- perr with n measurements
    T = studentT(n - 2)
    upper_tstat = (p - tolerance) / perr * np.sqrt(n)
    lower_tstat = (p + tolerance) / perr * np.sqrt(n)
    upper_p = (1 - T.cdf(abs(upper_tstat)))
    lower_p = (1 - T.cdf(abs(lower_tstat)))
    return lower_p, upper_p

x = np.arange(500)
yerr = np.ones_like(x) * 0.01
y = np.random.normal(0, yerr)
p, pcov = curve_fit(lambda r, m, c: c + (m*r), x, y, None, yerr, absolute_sigma=True)
perr = np.sqrt(np.diag(pcov))
pvalues = studentt_test(p[0], perr[0], 0.01, len(x))

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.