Sep-09-2020, 03:08 PM
Hi, I'm very new to pyhton and I've found the attached code from this link https://www.neuralengine.org/res/kernel.html#Code
I would like to run it on Anaconda3 Spyder 3.8.
The code does Kernel Bandwidth Optimization for a specific array. I've assigned x to some decimal numbers and run the function as below:
The whole code is below:
I would like to run it on Anaconda3 Spyder 3.8.
The code does Kernel Bandwidth Optimization for a specific array. I've assigned x to some decimal numbers and run the function as below:
x={0.30,0.30,0.31,0.32,0.30,0.5,0.6,0.8,0.8,0.81,0.82,0.80,0.81,1.0,1.3} ssvkernel(x,None,15,15,'Gauss')But it gives below errors.
Error:Traceback (most recent call last):
File "<ipython-input-8-4eb7a6b30c47>", line 1, in <module>
ssvkernel(x,None,15,15,'Gauss')
File "C:\Users\Ersin\Documents\ssvkernel.py", line 77, in ssvkernel
dx = np.sort(np.diff(np.sort(x)))
File "<__array_function__ internals>", line 5, in sort
File "C:\ProgramData\Anaconda3\lib\site-packages\numpy\core\fromnumeric.py", line 989, in sort
a.sort(axis=axis, kind=kind, order=order)
AxisError: axis -1 is out of bounds for array of dimension 0
Can some one help me to figure out the issue?The whole code is below:
import numpy as np def ssvkernel(x, tin=None, M=80, nbs=1e2, WinFunc='Boxcar'): """ Generates a locally adaptive kernel-density estimate for one-dimensional data. The user provides a one-dimensional vector of samples drawn from some underlying unknown distribution, and optionally the values where they want to estimate the probability density of that distribution. The algorithm solves an optimization problem to identify variable bandwidths across the domain where the data is provided. The optimization is based on a principle of minimizing expected L2 loss function between the kernel estimate and an unknown underlying density function. An assumption is merely that samples are drawn from the density independently of each other. The locally adaptive bandwidth is obtained by iteratively computing optimal fixed-size bandwidths wihtihn local intervals. The optimal bandwidths are selected such that they are selected in the intervals that are gamma times larger than the optimal bandwidths themselves. The paramter gamma is optimized by minimizing the L2 risk estimate. Parameters ---------- x : array_like The one-dimensional samples drawn from the underlying density tin : array_like, optional The values where the density estimate is to be evaluated in generating the output 'y'. Default value = None. M : int, optional The number of window sizes to evaluate. Default value = 80. nbs : int, optional The number of bootstrap samples to use in estimating the [0.05, 0.95] confidence interval of the output 'y'. WinFunc : string, optional The type of window function to use in estimating local bandwidth. Choose from one of 'Boxcar', 'Laplace', 'Cauchy' and 'Gauss'. Default value = 'Gauss'. Returns ------- y : array_like The estimated density, evaluated at points t / tin. t : array_like The points where the density estimate 'y' is evaluated. optw : array_like The optimal local kernel bandwidths at 't'. gs : array_like The stiffness constants of the variables bandwidths evaluated. C : array_like Cost functions associated with stiffness constraints. confb95 : array_like The 5% and 95% confidence interval of the kernel density estimate 'y'. Has dimensions 2 x len(y). confb95[0,:] corresponds to the 5% interval, and confb95[1,:] corresponds to the 95% interval. yb : array_like The bootstrap samples used in estimating confb95. Each row corresponds to one bootstrap sample. See Also -------- sshist, sskernel References ---------- .. [1] H. Shimazaki and S. Shinomoto, "Kernel Bandwidth Optimization in Spike Rate Estimation," in Journal of Computational Neuroscience 29(1-2): 171–182, 2010 http://dx.doi.org/10.1007/s10827-009-0180-4 """ # set argument 't' if not provided if tin is None: T = np.max(x) - np.min(x) dx = np.sort(np.diff(np.sort(x))) dt_samp = dx[np.nonzero(dx)][0] tin = np.linspace(np.min(x), np.max(x), min(np.ceil(T / dt_samp), 1e3)) t = tin x_ab = x[(x >= min(tin)) & (x <= max(tin))] else: T = np.max(x) - np.min(x) x_ab = x[(x >= min(tin)) & (x <= max(tin))] dx = np.sort(np.diff(np.sort(x))) dt_samp = dx[np.nonzero(dx)][0] if dt_samp > min(np.diff(tin)): t = np.linspace(min(tin), max(tin), min(np.ceil(T / dt_samp, 1e3))) else: t = tin # calculate delta t dt = min(np.diff(t)) # create the finest histogram thist = np.concatenate((t, (t[-1]+dt)[np.newaxis])) y_hist = np.histogram(x_ab, thist-dt/2)[0] / dt L = y_hist.size N = sum(y_hist * dt).astype(np.float) # initialize window sizes W = logexp(np.linspace(ilogexp(5 * dt), ilogexp(T), M)) # compute local cost functions c = np.zeros((M, L)) for j in range(M): w = W[j] yh = fftkernel(y_hist, w / dt) c[j, :] = yh**2 - 2 * yh * y_hist + 2 / (2 * np.pi)**0.5 / w * y_hist # initialize optimal ws optws = np.zeros((M, L)) for i in range(M): Win = W[i] C_local = np.zeros((M, L)) for j in range(M): C_local[j, :] = fftkernelWin(c[j, :], Win / dt, WinFunc) n = np.argmin(C_local, axis=0) optws[i, :] = W[n] # golden section search for stiffness parameter of variable bandwidths k = 0 gs = np.zeros((30, 1)) C = np.zeros((30, 1)) tol = 1e-5 a = 1e-12 b = 1 phi = (5**0.5 + 1) / 2 c1 = (phi - 1) * a + (2 - phi) * b c2 = (2 - phi) * a + (phi - 1) * b f1 = CostFunction(y_hist, N, t, dt, optws, W, WinFunc, c1)[0] f2 = CostFunction(y_hist, N, t, dt, optws, W, WinFunc, c2)[0] while (np.abs(b-a) > tol * (abs(c1) + abs(c2))) & (k < 30): if f1 < f2: b = c2 c2 = c1 c1 = (phi - 1) * a + (2 - phi) * b f2 = f1 f1, yv1, optwp1 = CostFunction(y_hist, N, t, dt, optws, W, WinFunc, c1) yopt = yv1 / np.sum(yv1 * dt) optw = optwp1 else: a = c1 c1 = c2 c2 = (2 - phi) * a + (phi - 1) * b f1 = f2 f2, yv2, optwp2 = CostFunction(y_hist, N, t, dt, optws, W, WinFunc, c2) yopt = yv2 / np.sum(yv2 * dt) optw = optwp2 # capture estimates and increment iteration counter gs[k] = c1 C[k] = f1 k = k + 1 # discard unused entries in gs, C gs = gs[0:k] C = C[0:k] # estimate confidence intervals by bootstrapping nbs = np.asarray(nbs) yb = np.zeros((nbs, tin.size)) for i in range(nbs): Nb = np.random.poisson(lam=N) idx = np.random.randint(0, N, Nb) xb = x_ab[idx] thist = np.concatenate((t, (t[-1]+dt)[np.newaxis])) y_histb = np.histogram(xb, thist - dt / 2)[0] idx = y_histb.nonzero() y_histb_nz = y_histb[idx] t_nz = t[idx] yb_buf = np.zeros((L, )) for k in range(L): yb_buf[k] = np.sum(y_histb_nz * Gauss(t[k] - t_nz, optw[k])) / Nb yb_buf = yb_buf / np.sum(yb_buf * dt) yb[i, :] = np.interp(tin, t, yb_buf) ybsort = np.sort(yb, axis=0) y95b = ybsort[np.floor(0.05 * nbs), :] y95u = ybsort[np.floor(0.95 * nbs), :] confb95 = np.concatenate((y95b[np.newaxis], y95u[np.newaxis]), axis=0) # return outputs y = np.interp(tin, t, yopt) optw = np.interp(tin, t, optw) t = tin return y, t, optw, gs, C, confb95, yb def CostFunction(y_hist, N, t, dt, optws, WIN, WinFunc, g): L = y_hist.size optwv = np.zeros((L, )) for k in range(L): gs = optws[:, k] / WIN if g > np.max(gs): optwv[k] = np.min(WIN) else: if g < min(gs): optwv[k] = np.max(WIN) else: idx = np.max(np.nonzero(gs >= g)) optwv[k] = g * WIN[idx] # Nadaraya-Watson kernel regression optwp = np.zeros((L, )) for k in range(L): if WinFunc == 'Boxcar': Z = Boxcar(t[k]-t, optwv / g) elif WinFunc == 'Laplace': Z = Laplace(t[k]-t, optwv / g) elif WinFunc == 'Cauchy': Z = Cauchy(t[k]-t, optwv / g) else: # WinFunc == 'Gauss' Z = Gauss(t[k]-t, optwv / g) optwp[k] = np.sum(optwv * Z) / np.sum(Z) # speed-optimized baloon estimator idx = y_hist.nonzero() y_hist_nz = y_hist[idx] t_nz = t[idx] yv = np.zeros((L, )) for k in range(L): yv[k] = np.sum(y_hist_nz * dt * Gauss(t[k]-t_nz, optwp[k])) yv = yv * N / np.sum(yv * dt) # cost function of estimated kernel cg = yv**2 - 2 * yv * y_hist + 2 / (2 * np.pi)**0.5 / optwp * y_hist Cg = np.sum(cg * dt) return Cg, yv, optwp def fftkernel(x, w): # forward padded transform L = x.size Lmax = L + 3 * w n = 2 ** np.ceil(np.log2(Lmax)) X = np.fft.fft(x, n.astype(np.int)) # generate kernel domain f = np.linspace(0, n-1, n) / n f = np.concatenate((-f[0: n / 2 + 1], f[1: n / 2 - 1 + 1][::-1])) # evaluate kernel K = np.exp(-0.5 * (w * 2 * np.pi * f) ** 2) # convolve and transform back from frequency domain y = np.real(np.fft.ifft(X * K, n)) y = y[0:L] return y def fftkernelWin(x, w, WinFunc): # forward padded transform L = x.size Lmax = L + 3 * w n = 2 ** np.ceil(np.log2(Lmax)) X = np.fft.fft(x, n.astype(np.int)) # generate kernel domain f = np.linspace(0, n-1, n) / n f = np.concatenate((-f[0: n / 2 + 1], f[1: n / 2 - 1 + 1][::-1])) t = 2 * np.pi * f # determine window function - evaluate kernel if WinFunc == 'Boxcar': a = 12**0.5 * w K = 2 * np.sin(a * t / 2) / (a * t) K[0] = 1 elif WinFunc == 'Laplace': K = 1 / (1 + (w * 2 * np.pi * f)**2 / 2) elif WinFunc == 'Cauchy': K = np.exp(-w * np.abs(2 * np.pi * f)) else: # WinFunc == 'Gauss' K = np.exp(-0.5 * (w * 2 * np.pi * f)**2) # convolve and transform back from frequency domain y = np.real(np.fft.ifft(X * K, n)) y = y[0:L] return y def Gauss(x, w): y = 1 / (2 * np.pi)**2 / w * np.exp(-x**2 / 2 / w**2) return y def Laplace(x, w): y = 1 / 2**0.5 / w * np.exp(-(2**0.5) / w / np.abs(x)) return y def Cauchy(x, w): y = 1 / (np.pi * w * (1 + (x / w)**2)) return y def Boxcar(x, w): a = 12**0.5 * w y = 1 / a y[np.abs(x) > a / 2] = 0 return y def logexp(x): y = np.zeros(x.shape) y[x < 1e2] = np.log(1+np.exp(x[x < 1e2])) y[x >= 1e2] = x[x >= 1e2] return y def ilogexp(x): y = np.zeros(x.shape) y[x < 1e2] = np.log(np.exp(x[x < 1e2]) - 1) y[x >= 1e2] = x[x >= 1e2] return y