alpha kernel is not starting at zero

[grg2rsr]

Contrary to the definition in the docs, the alpha kernel does not start with non-negative values only after t = 0 in elephant 0.4.1.

to demonstrate

import neo
import elephant as ele
import scipy as sp
import matplotlib.pyplot as plt
import quantities as pq

st = neo.SpikeTrain([1]*pq.s,t_start=0*pq.s,t_stop=2*pq.s)
kernel = ele.kernels.AlphaKernel(200*pq.ms)
fs = 0.1 * pq.ms
asig = ele.statistics.instantaneous_rate(st,t_start=st.t_start,t_stop=st.t_stop,sampling_period=fs,kernel=kernel)
plt.plot(asig.times.rescale(pq.s),asig)
[plt.axvline(t) for t in st.times]

produces

edit: which is almost the median (not exactly probably the kernel never decays back to zero). I looked at the AlphaKernel class in kernels.py and it looks fine, is there some median centering of the kernel going on somewhere else, that might need to be disabled?

[grg2rsr]

found in statistics.instantaneous_rate() line 757

    if not trim:
        r = r[kernel.median_index(t_arr):-(kernel(t_arr).size -
                                           kernel.median_index(t_arr))]

what is the reason for this line?

[mdenker]

This is an interesting use of instantaneous_rate with a non-symmetric kernel that I had not thought of. Indeed, this plot looks a bit strange at first. As you point out, the asymmetric kernel itself is not centered, but is 0 for x<0, and the centering is only performed by the instantaneous rate function. I assume the reason for centering the kernel on the mean is to have equal weight of the spike in the past and the future. This way the time of spike occurrence indicates the median of the probability mass. As such, I think that this choice makes sense as default behavior when interpreting the output of the function as an instantaneous rate.

In your example, I assume you would like to have the alpha kernel starting at the spike time, e.g., to obtain something like an EPSP after the spike?

In this case, one could either introduce a parameter like center=True to the instantaneous_rate, which will not center the kernel to the median when set to False. Or alternatively, to have a separate function convolve_spiketrain(), which does this?

[grg2rsr]

In your example, I assume you would like to have the alpha kernel starting at the spike time, e.g., to obtain something like an EPSP after the spike?

No not really, just normal firing rate estimation, but I chose an asymmetric kernel starting at t=0 to avoid a rise in firing rate before the response onset, so to say to avoid 'looking into the future' from the perspective of the neuron.
I think a center=True keyword makes sense to have, as the kernel itself was as I expected, but the output of the function was not. An extra function would not be necessary.

[mdenker]

Ah, ok, now I get it. I have created a branch with a potential fix at
https://github.com/INM-6/elephant/tree/feature/non-symmetric_instantaneous_rate

I think that this works correctly for all combinations of True/False for trim and the new option center_kernel using a modified version of your example:

import elephant as ele
import scipy as sp
import matplotlib.pyplot as plt
import quantities as pq
reload(elephant.statistics)
st = neo.SpikeTrain([3]*pq.s,t_start=0*pq.s,t_stop=10*pq.s)
kernel = ele.kernels.AlphaKernel(200*pq.ms)
fs = 0.1 * pq.ms
asig = ele.statistics.instantaneous_rate(st,t_start=st.t_start,t_stop=st.t_stop,sampling_period=fs,kernel=kernel,trim=True,center_kernel=True)
plt.plot(asig.times.rescale(pq.s),asig)
[plt.axvline(t) for t in st.times]
plt.show()

Could you please have a look if this addresses your issue? If so, I would PR this to the main branch (and write a short unit test for it).

[grg2rsr]

This looks fine. Thank you very much!

[grg2rsr]

I just stumbled over this - again. And then I found my bug report from a few years ago. As it seems, it was fixed, but now I get

TypeError: instantaneous_rate() got an unexpected keyword argument 'center_kernel'

was this never merged into the main branch, and the bugfix got lost?

[jpgill86]

I'd love to see the center_kernel feature added as well. I put together my own workaround for using the AlphaKernel to represent EPSPs sometime ago ("CausalAlphaKernel", code, docs), but it is very much a hack and far from ideal! (And had I noticed this issue and @mdenker's feature branch before, I probably would have just used it instead!)

[grg2rsr]

@jpgill86 very cool, thank you for posting this workaround

[dizcza]

I'll work on the issue in conjunction with #288.
Firstly, I need to understand the current status of kernels in Elephant and median_index() function in particular.
After this, I'll include the center_kernel parameter in kernels.
It'll take time though.

[dizcza]

center_kernel argument is now in the instantaneous_rate() function (True by default) in the master branch. Feel free to test it and give us feedback.