|
| 1 | +# -*- coding: utf-8 -*- |
| 2 | +""" |
| 3 | +Tools for analyzing the autocorrelation of a time series |
| 4 | +""" |
| 5 | +import numpy as np |
| 6 | + |
| 7 | +def autocorrfxn(timeseries,lagmax): |
| 8 | + ts = np.asarray(timeseries) |
| 9 | + ts -= np.average(ts) # Set to mean 0 |
| 10 | + N = len(timeseries) |
| 11 | + corrfxn = np.zeros(lagmax) |
| 12 | + for dt in xrange(lagmax): |
| 13 | + corrfxn[dt] = (np.dot(timeseries[0:N-dt],timeseries[dt:N])) # sum of ts[t+dt]*ts[t] |
| 14 | + corrfxn /= corrfxn[0] # Normalize |
| 15 | + return corrfxn |
| 16 | + |
| 17 | + |
| 18 | +def ipce(timeseries,lagmax=None): |
| 19 | + """ |
| 20 | + Initial positive correlation time estimator |
| 21 | + """ |
| 22 | + if lagmax == None: |
| 23 | + lagmax = len(timeseries)/2 |
| 24 | + corrfxn = autocorrfxn(timeseries,lagmax) |
| 25 | + i = 0 |
| 26 | + t = 0 |
| 27 | + while i < 0.5*lagmax: |
| 28 | + gamma = corrfxn[2*i] + corrfxn[2*i+1] |
| 29 | + if gamma < 0.0: |
| 30 | + print 'stop at ',2*i |
| 31 | + break |
| 32 | + else: |
| 33 | + t += gamma |
| 34 | + i += 1 |
| 35 | + tau = 2*t - 1 |
| 36 | + std = np.std(timeseries) |
| 37 | + mean = np.average(timeseries) |
| 38 | + sigma = std * tau / len(timeseries) |
| 39 | + return tau, mean, sigma |
| 40 | + |
| 41 | +def icce(timeseries,lagmax=None): |
| 42 | + """ |
| 43 | + Initial convex correlation time estimator |
| 44 | + REWRITE TO BE MORE EFFICIENT! |
| 45 | + """ |
| 46 | + if lagmax == None: |
| 47 | + lagmax = len(timeseries)/2 |
| 48 | + corrfxn = autocorrfxn(timeseries,lagmax) |
| 49 | + t = corrfxn[0] + corrfxn[1] |
| 50 | + i = 1 |
| 51 | + gammapast = t |
| 52 | + gamma = corrfxn[2*i] = corrfxn[2*i+1] |
| 53 | + while i < 0.5*lagmax: |
| 54 | + gammafuture = corrfxn[2*i+2] + corrfxn[2*i+3] |
| 55 | + if gamma > 0.5*(gammapast+gammafuture) : |
| 56 | + print 'stop at ',2*i |
| 57 | + break |
| 58 | + else: |
| 59 | + t += gamma |
| 60 | + gammapast = gamma |
| 61 | + gamma = gammafuture |
| 62 | + i += 1 |
| 63 | + tau = 2*t - 1 |
| 64 | + std = np.std(timeseries) |
| 65 | + mean = np.average(timeseries) |
| 66 | + sigma = std * tau / len(timeseries) |
| 67 | + return tau, mean, sigma |
| 68 | + |
0 commit comments