from numpy import * from matplotlib.pyplot import * # Cross Correlation to compute ACF def crosscorr(A, B): ftA = fft.rfft(A) ftB = fft.rfft(B) return 1./len(A)*fft.irfft(ftA.conj()*ftB) # Binning print "Binning Analysis" def block_mean(x, k): """Computes block averages of blocks of size k of the data series x.""" N = len(x) Nb = int(N/k) blocks = empty(Nb) for i in range(Nb): blocks[i] = x[i*k:(i+1)*k].mean() return blocks def binning(x, kmax): tau_bins = [] errs = [] sigma2o = x.var() for k in range(1, kmax+1): blocks = block_mean(x, k) Nb = len(blocks) sigma2b = blocks.var() tau_bins.append(0.5*sigma2b*k/sigma2o) errs.append(sqrt(sigma2b/Nb)) return array(tau_bins), array(errs) d0, d1, d2, d3, d4 = load('artificial.npy') figure('Artificial Data') subplot(221, title='Data') plot(d0[:1000]) plot(d1[:1000]) plot(d2[:1000]) plot(d3[:1000]) plot(d4[:1000]) # ACF d0 -= d0.mean() d1 -= d1.mean() d2 -= d2.mean() d3 -= d3.mean() d4 -= d4.mean() acf0 = crosscorr(d0, d0) acf1 = crosscorr(d1, d1) acf2 = crosscorr(d2, d2) acf3 = crosscorr(d3, d3) acf4 = crosscorr(d4, d4) subplot(222, title="ACF") plot(acf0[:1000]) plot(acf1[:1000]) plot(acf2[:1000]) plot(acf3[:1000]) plot(acf4[:1000]) # Binning print "Compute binning..." print " 0" tau0, errs0 = binning(d0, 1000) print " 1" tau1, errs1 = binning(d1, 1000) print " 2" tau2, errs2 = binning(d2, 1000) print " 3" tau3, errs3 = binning(d3, 1000) print " 4" tau4, errs4 = binning(d4, 1000) subplot(223, title='Binning tau') plot(tau0) plot(tau1) plot(tau2) plot(tau3) plot(tau4) subplot(224, title='Binning error') plot(errs0) plot(errs1) plot(errs2) plot(errs3) plot(errs4) tight_layout() # Simulation Data t, n0, n1, n2 = load('simulation.npy') t_eq = 20000 n_eq=0 while t[n_eq] < t_eq: n_eq += 1 figure('Simulation Data') subplot(221, title='Data') plot(n0[n_eq:]) plot(n1[n_eq:]) plot(n2[n_eq:]) # ACF d0 = n0[n_eq:]-n0[n_eq:].mean() d0 /= d0.std() acf0 = crosscorr(d0, d0) d1 = n1[n_eq:]-n1[n_eq:].mean() d1 /= d1.std() acf1 = crosscorr(d1, d1) d2 = n2[n_eq:]-n2[n_eq:].mean() d2 /= d2.std() acf2 = crosscorr(d2, d2) subplot(222, title="ACF") plot(acf0[:1000]) plot(acf1[:1000]) plot(acf2[:1000]) # Binning print "Compute binning..." print " 0" tau0, errs0 = binning(d0, 1000) print " 1" tau1, errs1 = binning(d1, 1000) print " 2" tau2, errs2 = binning(d2, 1000) subplot(223, title='Binning tau') plot(tau0) plot(tau1) plot(tau2) subplot(224, title='Binning error') plot(errs0) plot(errs1) plot(errs2) tight_layout() show()