Fixes NeuralEnsemble#245, where the automatic kernel selection yields…

… incorrect values (NeuralEnsemble#246)

* Fixes NeuralEnsemble#245, where the automatic kernel selection yields incorrect values

* Fixed spelling error, and remove irrelevant bootstrapping

* Also removed bootstrapping from regression test

* fixed python2 ssl CERTIFICATE_VERIFY_FAILED when downloading a file over https

* fixed instantaneous_rate numeric precision warning

elephant/statistics.py

@@ -535,152 +535,6 @@ def make_kernel(form, sigma, sampling_period, direction=1):  # pragma: no cover
     return kernel, norm, m_idx
 
 
-# to finally be taken out of Elephant
-def oldfct_instantaneous_rate(spiketrain, sampling_period, form,
-                              sigma='auto', t_start=None, t_stop=None,
-                              acausal=True, trim=False):  # pragma: no cover
-    """
-    Estimate instantaneous firing rate by kernel convolution.
-
-    Parameters
-    -----------
-    spiketrain: 'neo.SpikeTrain'
-        Neo object that contains spike times, the unit of the time stamps
-        and t_start and t_stop of the spike train.
-    sampling_period : Quantity
-        time stamp resolution of the spike times. the same resolution will
-        be assumed for the kernel
-    form : {'BOX', 'TRI', 'GAU', 'EPA', 'EXP', 'ALP'}
-        Kernel form. Currently implemented forms are BOX (boxcar),
-        TRI (triangle), GAU (gaussian), EPA (epanechnikov), EXP (exponential),
-        ALP (alpha function). EXP and ALP are asymmetric kernel forms and
-        assume optional parameter `direction`.
-    sigma : string or Quantity
-        Standard deviation of the distribution associated with kernel shape.
-        This parameter defines the time resolution of the kernel estimate
-        and makes different kernels comparable (cf. [1] for symmetric kernels).
-        This is used here as an alternative definition to the cut-off
-        frequency of the associated linear filter.
-        Default value is 'auto'. In this case, the optimized kernel width for
-        the rate estimation is calculated according to [1]. Note that the
-        automatized calculation of the kernel width ONLY works for gaussian
-        kernel shapes!
-    t_start : Quantity (Optional)
-        start time of the interval used to compute the firing rate, if None
-        assumed equal to spiketrain.t_start
-        Default:None
-    t_stop : Qunatity
-        End time of the interval used to compute the firing rate (included).
-        If none assumed equal to spiketrain.t_stop
-        Default:None
-    acausal : bool
-        if True, acausal filtering is used, i.e., the gravity center of the
-        filter function is aligned with the spike to convolve
-        Default:None
-    m_idx : int
-        index of the value in the kernel function vector that corresponds
-        to its gravity center. this parameter is not mandatory for
-        symmetrical kernels but it is required when asymmetrical kernels
-        are to be aligned at their gravity center with the event times if None
-        is assumed to be the median value of the kernel support
-        Default : None
-    trim : bool
-        if True, only the 'valid' region of the convolved
-        signal are returned, i.e., the points where there
-        isn't complete overlap between kernel and spike train
-        are discarded
-        NOTE: if True and an asymmetrical kernel is provided
-        the output will not be aligned with [t_start, t_stop]
-
-    Returns
-    -------
-    rate : neo.AnalogSignal
-        Contains the rate estimation in unit hertz (Hz).
-        Has a property 'rate.times' which contains the time axis of the rate
-        estimate. The unit of this property is the same as the resolution that
-        is given as an argument to the function.
-
-    Raises
-    ------
-    TypeError:
-        If argument value for the parameter `sigma` is not a quantity object
-        or string 'auto'.
-
-    See also
-    --------
-    elephant.statistics.make_kernel
-
-    References
-    ----------
-    ..[1] H. Shimazaki, S. Shinomoto, J Comput Neurosci (2010) 29:171–182.
-    """
-    warnings.simplefilter('always', DeprecationWarning)
-    warnings.warn("deprecated", DeprecationWarning, stacklevel=2)
-    if sigma == 'auto':
-        form = 'GAU'
-        unit = spiketrain.units
-        kernel_width = sskernel(spiketrain.magnitude, tin=None,
-                                bootstrap=True)['optw']
-        sigma = kw2sigma(form) * kernel_width * unit
-    elif not isinstance(sigma, pq.Quantity):
-        raise TypeError('sigma must be either a quantities object or "auto".'
-                        ' Found: %s, value %s' % (type(sigma), str(sigma)))
-
-    kernel, norm, m_idx = make_kernel(form=form, sigma=sigma,
-                                      sampling_period=sampling_period)
-    units = pq.CompoundUnit(
-        "%s*s" % str(sampling_period.rescale('s').magnitude))
-    spiketrain = spiketrain.rescale(units)
-    if t_start is None:
-        t_start = spiketrain.t_start
-    else:
-        t_start = t_start.rescale(spiketrain.units)
-
-    if t_stop is None:
-        t_stop = spiketrain.t_stop
-    else:
-        t_stop = t_stop.rescale(spiketrain.units)
-
-    time_vector = np.zeros(int((t_stop - t_start)) + 1)
-
-    spikes_slice = spiketrain.time_slice(t_start, t_stop) \
-        if len(spiketrain) else np.array([])
-
-    for spike in spikes_slice:
-        index = int((spike - t_start))
-        time_vector[index] += 1
-
-    r = norm * scipy.signal.fftconvolve(time_vector, kernel, 'full')
-    if np.any(r < 0):
-        warnings.warn('Instantaneous firing rate approximation contains '
-                      'negative values, possibly caused due to machine '
-                      'precision errors')
-
-    if acausal:
-        if not trim:
-            r = r[m_idx:-(kernel.size - m_idx)]
-
-        elif trim:
-            r = r[2 * m_idx:-2 * (kernel.size - m_idx)]
-            t_start = t_start + m_idx * spiketrain.units
-            t_stop = t_stop - ((kernel.size) - m_idx) * spiketrain.units
-
-    else:
-        if not trim:
-            r = r[m_idx:-(kernel.size - m_idx)]
-
-        elif trim:
-            r = r[2 * m_idx:-2 * (kernel.size - m_idx)]
-            t_start = t_start + m_idx * spiketrain.units
-            t_stop = t_stop - ((kernel.size) - m_idx) * spiketrain.units
-
-    rate = neo.AnalogSignal(signal=r.reshape(r.size, 1),
-                            sampling_period=sampling_period,
-                            units=pq.Hz, t_start=t_start)
-
-    return rate, sigma
-
-
 def instantaneous_rate(spiketrain, sampling_period, kernel='auto',
                        cutoff=5.0, t_start=None, t_stop=None, trim=False):
     """
@@ -702,13 +556,13 @@ def instantaneous_rate(spiketrain, sampling_period, kernel='auto',
         The kernel is used for convolution with the spike train and its
         standard deviation determines the time resolution of the instantaneous
         rate estimation.
-        Default: 'auto'. In this case, the optimized kernel width for the 
+        Default: 'auto'. In this case, the optimized kernel width for the
         rate estimation is calculated according to [1] and with this width
-        a gaussian kernel is constructed. Automatized calculation of the 
+        a gaussian kernel is constructed. Automatized calculation of the
         kernel width is not available for other than gaussian kernel shapes.
     cutoff : float
         This factor determines the cutoff of the probability distribution of
-        the kernel, i.e., the considered width of the kernel in terms of 
+        the kernel, i.e., the considered width of the kernel in terms of
         multiples of the standard deviation sigma.
         Default: 5.0
     t_start : Time Quantity (optional)
@@ -734,8 +588,9 @@ def instantaneous_rate(spiketrain, sampling_period, kernel='auto',
     Returns
     -------
     rate : neo.AnalogSignal
-        Contains the rate estimation in unit hertz (Hz).
-        Has a property 'rate.times' which contains the time axis of the rate
+        Contains the rate estimation in unit hertz (Hz). In case a list of
+        spike trains was given, this is the combined rate of all spike trains
+        (not the average rate). 'rate.times' contains the time axis of the rate
         estimate. The unit of this property is the same as the resolution that
         is given via the argument 'sampling_period' to the function.
 
@@ -794,18 +649,13 @@ def instantaneous_rate(spiketrain, sampling_period, kernel='auto',
         raise ValueError("The sampling period must be larger than zero.")
 
     if kernel == 'auto':
-        kernel_width = sskernel(spiketrain.magnitude, tin=None,
-                                bootstrap=True)['optw']
-        if kernel_width is None:
+        kernel_width_sigma = sskernel(
+            spiketrain.magnitude, tin=None, bootstrap=False)['optw']
+        if kernel_width_sigma is None:
             raise ValueError(
                 "Unable to calculate optimal kernel width for "
                 "instantaneous rate from input data.")
-        unit = spiketrain.units
-        sigma = 1 / (2.0 * 2.7) * kernel_width * unit
-        # factor 2.0 connects kernel width with its half width,
-        # factor 2.7 connects half width of Gaussian distribution with
-        #             99% probability mass with its standard deviation.
-        kernel = kernels.GaussianKernel(sigma)
+        kernel = kernels.GaussianKernel(kernel_width_sigma * spiketrain.units)
     elif not isinstance(kernel, kernels.Kernel):
         raise TypeError(
             "kernel must be either instance of :class:`Kernel` "
@@ -863,7 +713,7 @@ def instantaneous_rate(spiketrain, sampling_period, kernel='auto',
 
     r = scipy.signal.fftconvolve(time_vector,
                                  kernel(t_arr).rescale(pq.Hz).magnitude, 'full')
-    if np.any(r < 0):
+    if np.any(r < -1e-8):  # abs tolerance in np.isclose
         warnings.warn("Instantaneous firing rate approximation contains "
                       "negative values, possibly caused due to machine "
                       "precision errors.")
@@ -1134,14 +984,16 @@ def cost_function(x, N, w, dt):
 def sskernel(spiketimes, tin=None, w=None, bootstrap=False):
     """
 
-    Calculates optimal fixed kernel bandwidth.
+    Calculates optimal fixed kernel bandwidth, given as the standard deviation
+    sigma.
 
     spiketimes: sequence of spike times (sorted to be ascending).
 
-    tin: (optional) time points at which the kernel bandwidth is to be estimated.
+    tin: (optional) time points at which the kernel bandwidth is to be
+    estimated.
 
-    w: (optional) vector of kernel bandwidths. If specified, optimal
-    bandwidth is selected from this.
+    w: (optional) vector of kernel bandwidths (standard deviation sigma). If
+    specified, optimal bandwidth is selected from this.
 
     bootstrap (optional): whether to calculate the 95% confidence
     interval. (default False)
@@ -1152,8 +1004,8 @@ def sskernel(spiketimes, tin=None, w=None, bootstrap=False):
 
     'y': estimated density,
     't': points at which estimation was computed,
-    'optw': optimal kernel bandwidth,
-    'w': kernel bandwidths examined,
+    'optw': optimal kernel bandwidth given as standard deviation sigma
+    'w': kernel bandwidths examined (standard deviation sigma),
     'C': cost functions of w,
     'confb95': (lower bootstrap confidence level, upper bootstrap confidence level),
     'yb': bootstrap samples.