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 0The 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
