AllanTools 2016.11 - now with Confidence Intervals!

at2016-11_wtih_ci

I've added confidence interval estimation to allantools, based on a 2004 paper by Greenhall & Riley: "Uncertainty of stability variances based on finite differences"

So far not much is automated so you have to run everything manually. After a normal call to allantools.adev() to calculate the ADEV we loop through each (tau, adev) pair and first call allantools.edf_greenhall() to get the number of equivalent-degrees-of-freedom (EDF), and then evaluate a confidence interval with allantools.confidence_interval(). A knowledge or estimate of the noise type "alpha" is required for edf_greenhall() - here we just assume alpha=0.

This example is on github at: https://github.com/aewallin/allantools/blob/master/examples/ci_demo.py

import numpy
import matplotlib.pyplot as plt
 
import allantools as at
 
# this demonstrates how to calculate confidence intervals for ADEV
# using the algorithms from Greenhall2004
data_file = '../tests/phasedat/PHASE.DAT'
 
def read_datafile(filename):
    p = []
    with open(filename) as f:
        for line in f:
            if not line.startswith("#"):  # skip comments
                p.append(float(line))
    return p
 
# read input data from file
phase = read_datafile(data_file)
# normal ADEV computation, giving naive 1/sqrt(N) errors
(taus,devs,errs,ns) = at.adev(phase, taus='octave')
 
# Confidence-intervals for each (tau,adev) pair separately.
cis=[]
for (t,dev) in zip(taus,devs):
    # Greenhalls EDF (Equivalent Degrees of Freedom)
    # alpha     +2,...,-4   noise type, either estimated or known
    # d         1 first-difference variance, 2 allan variance, 3 hadamard variance
    #           we require: alpha+2*d >1     (is this ever false?)
    # m         tau/tau0 averaging factor
    # N         number of phase observations
    edf = at.edf_greenhall( alpha=0, d=2, m=t, N=len(phase), overlapping = False, modified=False )
    # with the known EDF we get CIs 
    # for 1-sigma confidence we set
    # ci = scipy.special.erf(1/math.sqrt(2)) = 0.68268949213708585
    (lo,hi) = at.confidence_intervals( dev=dev, ci=0.68268949213708585, edf=edf )
    cis.append( (lo,hi) )
 
# now we are ready to print and plot the results
print "Tau\tmin Dev\t\tDev\t\tMax Dev"
for (tau,dev,ci) in zip(taus,devs,cis):
    print "%d\t%f\t%f\t%f" % (tau, ci[0], dev, ci[1] )
""" output is
Tau	min Dev		Dev		Max Dev
1	0.285114	0.292232	0.299910
2	0.197831	0.205102	0.213237
4	0.141970	0.149427	0.158198
8	0.102541	0.110135	0.119711
16	0.056510	0.062381	0.070569
32	0.049153	0.056233	0.067632
64	0.027109	0.032550	0.043536
128	0.026481	0.033855	0.055737
256	0.007838	0.010799	0.031075
"""
plt.figure(figsize=(12,8))
plt.gca().set_yscale('log')
plt.gca().set_xscale('log')
err_lo = [ d-ci[0] for (d,ci) in zip(devs,cis)]
err_hi = [ ci[1]-d for (d,ci) in zip(devs,cis)]
 
plt.errorbar(taus, devs, yerr=[ err_lo, err_hi ] ,fmt='o')
plt.grid()
plt.xlabel('Tau (s)')
plt.ylabel('ADEV')
plt.title('AllanTools 2016.11 - now with Confidence Intervals!')
# just to check plot the intervals as dots also
plt.plot(taus, [ci[0] for ci in cis],'r.')
plt.plot(taus, [ci[1] for ci in cis],'g.')
 
plt.show()

scipy.stats.chi2.ppf() without scipy

Allantools has a dependence on scipy because it uses scipy.stats.chi2.ppf() (inverse of the chi-squared cumulative distribution function) in the code for confidence intervals.

When testing scipy takes quite a long time to install and the whole package seems a bit overkill for just this one function.

So I tried implementing it using just numpy. It sort of works, but not for corner cases where p is very close to 1.0 or k is a big number. I get for example:

    top = pow(x, s)*math.exp(-x)*pow(x,k)
OverflowError: (34, 'Numerical result out of range')

Maybe there's something trick when numerically evaluating lower_gamma() using the series expansion that I am missing...

import scipy.stats
import math
import numpy
 
def bisection(function, k, p, a, b, tol):
	# http://code.activestate.com/recipes/578417-bisection-method-in-python/
	assert (function(a,k)-p)*(function(b,k)-p) < 0 # a, b must bracket root
	c = (a+b)/2.0
	while (b-a)/2.0 > tol:
		if (function(c, k)-p) == 0:
			return c
		elif (function(a,k)-p)*(function(c,k)-p) < 0:
			b = c
		else :
			a = c
		c = (a+b)/2.0
	return c
 
def lower_gamma(s,x):
	# lower incomplete gamma function
	# https://en.wikipedia.org/wiki/Incomplete_gamma_function#Evaluation_formulae
	g  = 0
	last_g = 1.0
	done = False
	tol = 1.0e-6
	k=0
	while not done:
		top = pow(x, s)*math.exp(-x)*pow(x,k)
		bot = numpy.prod( [float(s+j) for j in range(k+1) ] )
		dg = float(top)/float(bot)
		if dg == float("Inf"):
			break
		g += dg
		k += 1
		if k>100: # get at least 100 terms in the sum
			if g==0:
				break
			delta = abs(dg/g)
			if delta == float("Inf"):
				break
			if delta < tol:
				done = True
		last_g = g
	return g
 
def chi2_cdf(x, k):
	# chi-squared cumulative density function
	# cdf(x; k) = lower_gamma(k/2, x/2) / gamma(k/2)
	return lower_gamma(k/2.0, x/2.0) / math.gamma(k/2.0)
 
def chi2_ppf(p, k):
	# chi-squared Percent point function (inverse of cdf percentiles).
	# look for x such that
	# p = chi2_cdf( x=chi2_ppf(p, k), k)
	tol = 1e-8
	lolim = 0
	hilim = k
	while (chi2_cdf(lolim,k)-p)*(chi2_cdf(hilim,k)-p) > 0: 
		hilim *= 1.5
	return bisection( chi2_cdf, k, p, lolim, hilim, tol)
 
print "scipy cdf: ",scipy.stats.chi2.cdf(55, 33)	
print "own   cdf: ",chi2_cdf(55, 33)
 
 
print "scipy ppf ", scipy.stats.chi2.ppf(0.4, 33)
print "  own ppf ", chi2_ppf(0.4, 33)
 
# test that we really found the inverse
print scipy.stats.chi2.cdf(scipy.stats.chi2.ppf(0.4, 33), 33)
print chi2_cdf( chi2_ppf(0.4, 33), 33 )
 
# try to check the scipy function against our own function
# for some random input of (p, k)
for n in range(100):
	k = numpy.random.randint(20, 200)
	p = numpy.random.random()
	print k, p,
	a=scipy.stats.chi2.ppf(p, k)
	b=chi2_ppf(p, k)
	ok = numpy.isclose(a, b)
	if ok:
		print  ok
	else:
		print  ok, a, b
	assert ok

Five colours of noise

Update: now with the colours matching in all graphs:

colorednoise

Time-series generated with colorednoise (following Kasdin&Walter), power-spectral-densities and Allan deviations computed with allantools, and compared to theoretical predictions in IEEE1139-2008.

colorednoise

The PSD lines and MDEV lines seem spot-on, but are the ADEV lines systematically a bit low?

Code here: example_noise_slopes.py

Hadamard total deviation in allantools

Following on from modified total deviation mtotdev() the Hadamard total deviation htotdev() algorithm is very similar but instead of phase data takes frequency data. Now included in allantools.

Here's a comarison against Stable32 and the NIST SP 1065 table values.

There's no bias-correction for now, and apparently it is customary not to use htotdev(m=1) at an averaging-factor of 1, but instead use ohdev(m=1). This is why the 'raw' value from allantools at tau=1 in the plot below is a factor of 2 too low.

htotdev_2016-03-27Both mtotdev() and htotdev() are very slow algorithms - help with making them faster would be appreciated!

Modified Total Deviation in allantools

I wrote a first rough (and very slow!) implementation of Modified Total Deviation mtotdev() for allantools.

mtotdev() combines the features of modified Allan deviation mdev() (being able to distinguish between white and flicker phase modulation) and Total deviation totdev() (better confidence intervals at large tau).

Here are some results. The Time Total Deviation ttotdev() follows trivially from this work also, since it is only a scaled version of mtotdev(). I used the "NBS14" 1000-point frequency dataset and compared my results against Stable32 and those in NIST SP 1065 (Table 31, page 108).

mtotdev_comparison_AW2016-03-24

ttotdev_comparison_AW2016-03-24

The figures show mtotdev() and ttotdev() from Stable32 runs and allantools. I've also added mdev() and tdev() traces to compare against. At first sight this looks strange, but mtotdev() is a biased estimator, and Stable32 applies a bias correction which explains the results.

Somewhat surprisingly Stable32 applies a bias-correction only when run in the "all-tau" mode. These numbers agree with those from the NIST SP 1065 table. For this dataset Stable32 is undecided on what power-law the dataset follows at large tau, which results in deviations that jump up and down (because a different bias-correction is applied at different tau, red datapoints).

When Stable32 is run in "octave-tau" or "decade-tau" mode no bias correction is applied. These numbers agree with the ones from allantools.mtotdev().

mtotdev_bias_AW2016-03-24

This figure shows the ratio of variances between allantools.mtotdev() (no bias correction) and a Stable32 all-tau run (bias correction applied). The lines correspond to the bias-correction values for different noise processes.

There also seems to be a misprint in NIST SP 1065 equation (28) page 26, where it says ttotdev() is scaled by tau-cubed (when tau-squared is correct). Who reads these things anyway - not many it seems 😉

Keysight 53230A noise floor test

We got a new 53230A counter to the lab, so I decided to run some basic tests on it.

I collected time interval data using a 1-PPS source (H-maser through a SRS DG645), and wired this with a T-connector from CH1 to CH2 with a ~1 m (10 ns delay) cable. This should show the noise floor for time interval measurements as well as CH1/CH2 timing skew when measured the other way around (i.e. from CH2 to CH1). The 10 MHz external reference (at the back) was connected to a H-maser.

53230A_PPS_skew
The results show standard deviations of 12 ps (CH1->CH2) and 11 ps (CH2->CH1) respecively, with a channel skew of 112 ps. Compare to the single-shot spec of sqrt(2)*20 ps = 28 ps and Agilent/Keysight's marketing video on youtube.

I also collected 10 MHz frequency counter readings on CH1 (source: H-maser) with gate times of 0.1 s, 1.0 s, and 10.0 s. I collected the data with a simple program that just calls the "READ?" function repeatedly, which does result in some dead-time between measurements.

Here are the results in terms of Allan deviation. I used allantools.

Keysight_53230A_noise_floor

The time interval noise floor looks like white phase noise with an Allan deviation of 1.8e-11/tau. This is consistent with the 12 ps RMS value found above. It is left as an exercise for the reader to show that ADEV(1s) = sqrt(3)*RMS-time-interval-noise (correct??).

The frequency counting noise floor depends on the gate time, and I get 5e-12/sqrt(tau), 2e-12/sqrt(tau), and 6e-13/sqrt(tau) for gate times of 0.1 s, 1.0 s, and 10.0 s, respectively. This looks like white frequency noise. Enrico Rubiola has notes on frequency counters that may explain the numbers.

White Rabbit Switch PPS-output test

We got some White Rabbit Switches and I did an initial test of the pulse-per-second (PPS) output stability. In contrast to earlier measurements that showed 200ps or so of white phase noise, the PPS output on the WRS now seems a lot more stable. For various reasons the noise-floor (red data) of our 53230A time-interval-counter is at around 50e-12 @ 1s, and the WRS PPS output is at very much the same level of stability. Another 53230A counter shows about 13 ps standard-deviation for a cable-delay measurement - so I may redo these measurements with that counter. In any case a real evaluation of the short-term stability requires a DMTD measurement at 10 MHz.