Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature/multifiltertest #142

Merged
merged 35 commits into from
Apr 3, 2018
Merged

Conversation

EmanueleLucrezia
Copy link
Contributor

This pull request adds the multiple filter test for the detection of rate change points in single unit spike trains. The algorithm is following the paper:

Messer, M., Kirchner, M., Schiemann, J., Roeper, J., Neininger, R., & Schneider, G. (2014). A multiple filter test for the detection of rate changes in renewal
processes with varying variance. The Annals of Applied Statistics, 8(4),2027-2067.

and the code is adapted from the R implementation:

Adapted from the published R implementation:
DOI: 10.1214/14-AOAS782SUPP;.r

@mdenker mdenker added the enhancement Editing an existing module, improving something label Mar 14, 2018
@mdenker
Copy link
Member

mdenker commented Mar 14, 2018

It seems the failing unit test -- which is related to reproducing published results -- is mainly failing as it depends on realization of the spike trains which are generated (with fixed CPs at 150, 180, and 500s). Although the random seed is fixed, the spike trains seems to differ on travis vs. when running it locally. In other words, the test data operates in a borderline regime, where the method will not always correctly detect all CPs in the data.

Thus, it may be a good idea to save the generated spike data explicitly and load it in the unit test.

In addition, one could remodel the current unit test such, that the differences between the rates are more pronounced such that the method is very sure to detect them.

@mdenker
Copy link
Member

mdenker commented Mar 14, 2018

In terms of file names -- would it be better to name the module change_point_detection? I think it's a more telling name, and would allow similar methods to be bundled under this name.

Independent of the module name decided upon, please name the corresponding unit test file test_XXX.py, where XXX is the module name.

@Junji110
Copy link
Contributor

Junji110 commented Mar 14, 2018

The unit test script (test_change_point_detection.py) imports elephant.multiple_filter_test and this raises an error. Please correct this line to the new module name.
Also, the example code in the docstring of change_point_detection.py imports multiple_filter_test. This also needs to be corrected.

@EmanueleLucrezia
Copy link
Contributor Author

The limit processes are not properly generated, namely, from pag. 10 of the article:
'Note that all limit processes (L h,t ) t are derived from the same Brownian motion in order to ensure comparability..'
I generated a different Brownian motions for every window h, instead to use always the same. This leads to a slightly higher threshold Q,

…brownian motion outside the for loop across windows
…en changed from 8, 13, 18, 16.5 Hz to 4, 13, 36, 16.5 Hz and the seed is not fixed anymore
@coveralls
Copy link
Collaborator

coveralls commented Mar 15, 2018

Coverage Status

Coverage increased (+1.2%) to 86.996% when pulling ea175dd on INM-6:feature/multifiltertest into 08d6ff1 on NeuralEnsemble:master.


>>> import quantities as pq
>>> import neo
>>> from elephant.multiple_filter_test import multiple_filter_test
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the module name (elephant.multiple_filter_test) needs to be corrected

variances of the limit process correspodning to `h`. This will be
used to normalize the `filter_process` in order to give to the every
maximum the same impact on the global statistic.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Explanation for test_quantile is missing, and the order of the arguments is not consistent with that of the definition of the function.

Returns:
--------
cps : list of list
one list for each `h`, containing the points detected with the
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is 'h'?

#print("detected point {0}".format(cp), "with filter {0}".format(h))
# before to repet the procedure the h-neighbourg of 'cp' detected
# is cut, because rate changes within it are explained by this cp
differences[np.where(
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This line is too complicated. Can be rewritten as:

            mask_fore = time_index > cp_index - int((h / dt_temp).simplified)
            mask_back = time_index < cp_index + int((h / dt_temp).simplified)
            differences[mask_fore & mask_back] = 0

# check if the neighbourhood of detected cp does not contain cps
# detected with other windows
neighbourhood_free = True
if i == 0:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This if-case is not necessary.
When i == 0, the for-loop in the else case actually doesn't iterate anything.

except ValueError:
raise ValueError("dt must be a time quantity")
t_fin_m = t_fin_sec.magnitude

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

t_in_sec and t_fin_sec are not actually used.
You should rather define them as t_in_sec = t_in.rescale(u).magnitude and t_fin_sec = t_fin.rescale(u).magnitude, and use them instead of t_in_m and t_fin_m.
This kind of confusing usage of rescale is seen in other functions as well.
Please correct them as well so that the treatment of units is consistent within the code.

matrix_normalized.append((matrix[i] - mean) / np.sqrt(var))
null_mean.append(mean)
null_var.append(var)
matrix_normalized = np.asanyarray(matrix_normalized)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This part of code is unnecessarily complicated and redundant.
The same result can be obtained by

null_mean = maxima_matrix.mean(axis=0)
null_var = maxima_matrix.var(axis=0)
matrix_normalized = (maxima_matrix - null_mean) / np.sqrt(null_var)

@alperyeg
Copy link
Contributor

@EmanueleLucrezia In the docstrings I inserted s at the headers, e.g Parameter -> Parameters. Otherwise, the documentation may not show these parts.

@mdenker
Copy link
Member

mdenker commented Mar 26, 2018

Comments by Messer:

  • It is particularly important to derive the different limit processes from the same underlying Brownian motion. The idea is that there is a single spike train on the 'side of real data’ that results in one single Brownian motion on the side of the ‘limit theory’. Then all Limit processes L are derived from the same Brownian motion as this corresponds to all filtered derivative processes G that are derived from the same spike train. (If you had different Brownian motions for different windows this would mean something like observing different spike trains in the same analysis, which does not make sense..)

  • The example in the original R-code is indeed chosen in a ‘borderline regime’ as you call it. The idea of that scenario was to show the necessity of different windows, i.e., besides large rate changes (small windows sufficient) choosing also a small rate increase where particularly a large window is necessary. Indeed, this results in the problem that the smaller rate increase is sometimes not detected (by chance). But of course feel free to increase the rates changes in you examples to ensure ‘robustness’ of their detection.

  • I also quickly scrolled along your code and by chance saw that you called $\alpha$ integer valued. I guess it should be floating as it is theoretically thought in the invertval (0,1)

@mdenker mdenker added this to the Release 0.5 milestone Mar 27, 2018
@alperyeg
Copy link
Contributor

@mdenker, see below our input to Messer's comments.

It is particularly important to derive the different limit processes from the same underlying Brownian motion (...)

This is addressed by commit 882d9ba

(...) alpha integer valued. I guess it should be floating as it is theoretically thought in the invertval (0,1)

With my latest commit alpha is now a floating number, also the documentation is accordingly updated.

(...) But of course feel free to increase the rates changes in you examples to ensure ‘robustness’ of their detection.

This is done with commit fd4b2da

@mdenker mdenker merged commit 1595c70 into NeuralEnsemble:master Apr 3, 2018
@Moritz-Alexander-Kern Moritz-Alexander-Kern deleted the feature/multifiltertest branch March 31, 2022 08:25
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement Editing an existing module, improving something
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants