master
Per A Brodtkorb 8 years ago
parent 07bcbde53b
commit 53a65e4f3d

@ -143,7 +143,8 @@ class CovarianceEstimator(object):
nfft = 2 ** nextpow2(n) nfft = 2 ** nextpow2(n)
raw_periodogram = abs(fft(x, nfft)) ** 2 / Ncens raw_periodogram = abs(fft(x, nfft)) ** 2 / Ncens
auto_cov = np.real(fft(raw_periodogram)) / nfft # ifft = fft/nfft since raw_periodogram is real! # ifft = fft/nfft since raw_periodogram is real!
auto_cov = np.real(fft(raw_periodogram)) / nfft
if self.flag.startswith('unbiased'): if self.flag.startswith('unbiased'):
# unbiased result, i.e. divide by n-abs(lag) # unbiased result, i.e. divide by n-abs(lag)

@ -862,7 +862,7 @@ class CyclePairs(PlotData):
# therefore we choose a slightly smaller bandwidth # therefore we choose a slightly smaller bandwidth
if aut_h == 1: if aut_h == 1:
h_norm = smoothcmat_hnorm(F,NOsubzero); h_norm = smoothcmat_hnorm(F, NOsubzero)
h = 0.7 * h_norm # Don't oversmooth h = 0.7 * h_norm # Don't oversmooth
# h0 = N^(-1/(d+4)); # h0 = N^(-1/(d+4));
@ -879,37 +879,37 @@ class CyclePairs(PlotData):
for i in range(n - 1): for i in range(n - 1):
for j in range(i + 1, n): for j in range(i + 1, n):
if F(i,j) != 0: if F[i, j] != 0:
F1 = exp(-1/(2*h**2)*((I-i)**2+(J-j)**2)); # Gaussian kernel F1 = exp(-1 / (2 * h**2) * ((I - i)**2 + (J - j)**2)) # Gaussian kernel
F1 = F1+F1.T; # Mirror kernel in diagonal F1 = F1 + F1.T # Mirror kernel in diagonal
F1 = np.triu(F1,1+NOsubzero); # Set to zero below and on diagonal F1 = np.triu(F1, 1 + NOsubzero) # Set to zero below and on diagonal
F1 = F[i, j] * F1 / np.sum(F1) # Normalize F1 = F[i, j] * F1 / np.sum(F1) # Normalize
Fsmooth = Fsmooth+F1; Fsmooth = Fsmooth + F1
# endif # endif
# endfor # endfor
# endfor # endfor
# endif method 1 or 2 # endif method 1 or 2
if method == 2: if method == 2:
Fpilot = Fsmooth/N; Fpilot = Fsmooth / N
Fsmooth = np.zeros(n,n); Fsmooth = np.zeros(n, n)
[I1,I2] = find(F>0); [I1, I2] = find(F > 0)
logg = 0; logg = 0
for i in range(len(I1)): # =1:length(I1): for i in range(len(I1)): # =1:length(I1):
logg = logg + F(I1(i),I2(i)) * log(Fpilot(I1(i),I2(i))); logg = logg + F(I1[i], I2[i]) * log(Fpilot(I1[i], I2[i]))
# endfor # endfor
g = np.exp(logg/N); g = np.exp(logg / N)
lamda = (Fpilot/g)**(-alpha); _lamda = (Fpilot / g)**(-alpha)
for i in range(n - 1): # = 1:n-1 for i in range(n - 1): # = 1:n-1
for j in range(i + 1, n): # = i+1:n for j in range(i + 1, n): # = i+1:n
if F[i, j] != 0: if F[i, j] != 0:
hi = h*lamda[i,j] hi = h * _lamda[i, j]
F1 = np.exp(-1/(2*hi**2)*((I-i)**2+(J-j)**2)); # Gaussian kernel F1 = np.exp(-1 / (2 * hi**2) * ((I - i)**2 + (J - j)**2)) # Gaussian kernel
F1 = F1+F1.T; # Mirror kernel in diagonal F1 = F1 + F1.T # Mirror kernel in diagonal
F1 = np.triu(F1,1+NOsubzero); # Set to zero below and on diagonal F1 = np.triu(F1, 1 + NOsubzero) # Set to zero below and on diagonal
F1 = F[i,j] * F1/np.sum(F1); # Normalize F1 = F[i, j] * F1 / np.sum(F1) # Normalize
Fsmooth = Fsmooth+F1; Fsmooth = Fsmooth + F1
# endif # endif
# endfor # endfor
# endfor # endfor
@ -1052,7 +1052,6 @@ class CyclePairs(PlotData):
else: else:
raise ValueError('Undefined discretization definition, ddef = {}'.format(ddef)) raise ValueError('Undefined discretization definition, ddef = {}'.format(ddef))
if np.any(ix): if np.any(ix):
cp1[ix], cp2[ix] = cp2[ix], cp1[ix] cp1[ix], cp2[ix] = cp2[ix], cp1[ix]
return np.asarray(cp1, type=int), np.asarray(cp2, type=int) return np.asarray(cp1, type=int), np.asarray(cp2, type=int)

Loading…
Cancel
Save