WIP: add reduction keyword to pwelch - #4000
Conversation
|
Sounds good.
Yes, private please.
That's what I would expect
You can just say that the output dimensions will be |
I would say just change an existing example to use |
Codecov Report
@@ Coverage Diff @@
## master #4000 +/- ##
==========================================
+ Coverage 83.01% 83.01% +<.01%
==========================================
Files 346 346
Lines 63074 63170 +96
Branches 9694 9715 +21
==========================================
+ Hits 52358 52439 +81
- Misses 8019 8028 +9
- Partials 2697 2703 +6 |
|
I've changed |
|
Why not make all well with the universe and do: |
|
Otherwise we'd need to add |
|
Scipy 0.12 problem: |
|
Did the old one operate on an axis or raveled array?
|
|
The old |
|
What do you suggest then? Backport to fixes?
|
|
Yes, I'd add this to fixes. I'll commit that later (tomorrow - European time) |
5d71401 to
2c6ad10
Compare
|
@Eric89GXL sorry for the delay, I've added |
|
One |
|
|
|
Ok, I'll rebase later today. |
75cebc6 to
f4483d6
Compare
|
@mmagnuski please set to MRG if it's ready from your end |
| scale = 1.0 / (fs * (win * win).sum()) | ||
| elif scaling == 'spectrum': | ||
| scale = 1.0 / win.sum()**2 | ||
| scale = 1.0 / win.sum() ** 2 |
There was a problem hiding this comment.
actually it's better to leave these so copy-paste diffs from scipy match
There was a problem hiding this comment.
ok, I'll revert that commit :)
| reduction is performed and psd's for individual windows are returned | ||
| (so that the output psds are of shape (windows, channels, frequencies) | ||
| for 2d input data and (windows, epochs, channels, frequencies) for 3d | ||
| input data) |
There was a problem hiding this comment.
Can just say "with the windows added as the first dimension"
There was a problem hiding this comment.
Oh, that's actually wrong. You previously commented that it should be -2 dim and I changed it that way. It makes sense as -2 - especially for epochs epochs x windows x channels x freqs)
| psds = np.concatenate(f_spectrogram, axis=0) | ||
| if reduction is not None: | ||
| reduction = _get_reduction(reduction) | ||
| psds = reduction(psds).reshape(np.hstack([dshape, -1])) |
There was a problem hiding this comment.
easier to use dshape + (-1,) than hstack (shapes should always be tuples)
| raise ValueError('reduction, if float, means proportion to trim in' | ||
| ' trimmed mean, which has to be > 0 and < 0.5, ' | ||
| 'got {}'.format(reduction)) | ||
| from .fixes import get_trim_mean |
There was a problem hiding this comment.
do you really need to nest this? fixes shouldn't require anything from utils...
|
Which example should I modify? I started changing the tutorial on freq analyses in channel space but it uses |
f4483d6 to
6b985ba
Compare
|
Travis errors seem to be unrelated: one job threw lots of docstyle errors from viz/utils at me and the other timed out when downloading mne test data. |
| else: | ||
| raise ValueError('reduction, if string, must be "mean" or "median"' | ||
| ', got {}'.format(reduction)) | ||
| elif isinstance(reduction, (float, np.float64, np.float32)): |
There was a problem hiding this comment.
import numbers
isinstance(reduction, numbers.Real)
There was a problem hiding this comment.
I still have to add this.
|
@mmagnuski yes I'll update plot_sensors_time_frequency.py with a section on the PSD of Raw data before going to Epochs. |
|
Sorry for the delays here, I was busy lately, I'll try to finish it this week. |
|
Will you use |
|
arfff my bad. Forget it
|
|
|
||
| # normalize multitaper and welch to put them on the same scale | ||
| psd_m /= psd_m.sum() / len(freqs_m) | ||
| psd_w /= psd_w.sum() / len(freqs_w) |
There was a problem hiding this comment.
psd_m.sum() -> psd_m.mean()
and why do you need to do this? both functions don't use the same spectral estimation?
There was a problem hiding this comment.
Maybe we don't - I was wondering whether this makes sense to add, but I thought that if almost everything in the tutorial was done with multitapering and now I am introducting Welch, a comparison would be nice. But I don't thinks it is necessary.
There was a problem hiding this comment.
please don't. Our PSD functions should be consistent in their normalizations
There was a problem hiding this comment.
Ok, now I see I didn't get your first question ("and why do you need to do this? both functions don't use the same spectral estimation?") - I thought you were asking why we should compare welch and multitapering in one figure. So yes, IIRC, there were big differences in psd returned by both functions - that's why I decided to normalize. I'll check that later to show an example of this.
| picks=[ch_index], fmin=2, fmax=17, n_jobs=1) | ||
|
|
||
| # drop channel dimension and average epochs | ||
| psd_m = psds_m.squeeze(axis=1).mean(axis=0) |
There was a problem hiding this comment.
please avoid squeeze. It's dangerous.
There was a problem hiding this comment.
Ok, I can use [:, 0] instead - but can squeeze be dangerous if I specify the axis?
There was a problem hiding this comment.
I prefer the [:, 0] it's more explicit
|
|
||
| ############################################################################### | ||
| # While in multitaper the averaging is done across independent realizations of | ||
| # the signal (using slepian tapers), welch method averages across time |
There was a problem hiding this comment.
slepian -> Slepian
welch -> Welch
|
|
||
| psd_w = psds_w[:, 0].mean(axis=0) | ||
| psd_w_trim = psds_w_trimmed[:, 0].mean(axis=0) | ||
| psd_w_med = psds_w_median[:, 0].mean(axis=0) |
There was a problem hiding this comment.
why do you need to average after a median?
what are you averaging?
There was a problem hiding this comment.
the epochs are averaged over
There was a problem hiding this comment.
don't you think it's weird to take the mean over epochs and the median of short windows within epochs?
There was a problem hiding this comment.
Hm, yes I agree that it is kind of weird. :) I can change it to median or I'll change that example later so that psd_welch is presented on raw data (as you suggested in other comments).
| # The reduction in power that can be seen in the figure above is due to the | ||
| # fact that values for power spectral density follow a positive skewed | ||
| # gamma-like distribution. Lets take a look at this distribution. First we will | ||
| # use ``combine=None`` to get all the welch windows without averaging. Notice |
There was a problem hiding this comment.
welch -> Welch
please check everywhere
|
see how it looks like:
https://4927-1301584-gh.circle-artifacts.com/0/tmp/circle-artifacts.JGVLnXY/html/auto_tutorials/plot_sensors_time_frequency.html
why would you want to do trimmed mean on EEG in the first place?
|
Well, I'm not sure if I would be doing that but in the original issue about adding reduction options to psd @kingjr proposed trimmed mean and your response was:
😋 |
|
(BTW - I checked the example before - it looks ok for me, is something wrong with how it looks?) |
|
the answer I was expecting is "I think it's a useful option to get robust
spectral estimates from potentially noisy data"
if you had said this, I would have asked that you illustrate this in an
example :)
|
|
(BTW - I checked the example before - it looks ok for me, is something
wrong with how it looks?)
it was just for others to see / comment
|
Ah, ok - I though you are asking why would I use trimmed mean over median.
Ok, I wasn't sure because I posted the link previously and the example didn't change since then. :) |
|
no rush :)
|
|
Hm, yes I agree that it is kind of weird. :) I can change it to median or
I'll change that example later so that psd_welch is presented on raw data
(as you suggested in other comments).
+1
just add a note here to crossref the two functions and mention that
psd_welch is available but is more natural to use on Raw cf
plot_visualize_raw.py
thanks a lot
|
|
@mmagnuski you need to rebase. any hope to complete this soon? I am thinking about pushing a release soon. |
|
Yes, I can take a look during the weekend and finish the remaining issues here. |
|
thanks great
|
|
@mmagnuski you need to rebase. Let's try to finish this before the release late august |
|
Yes, I'll finish this within a couple of days, sorry. |
|
Ping @mmagnuski it's getting close to release time, any time to polish this one? |
|
Yes, I should finally have time for that soon. Another option is that I revert changes to tutorial, we merge this now and I add a follow-up PR for changes in the tutorial. Everything else is finished here. |
|
revert changes in the tutorial you touched and just show the new option in:
https://martinos.org/mne/stable/auto_tutorials/plot_visualize_raw.html
we do raw.plot_psd there. Again let's show case this option for raw object
|
|
@mmagnuski any time for this one now? |
|
Yes, I should finally have time for this in next few days! :) |
|
@mmagnuski any time to come back to this one? You'll probably want to start fresh because of how many viz / PSD functions have changed :( |
larsoner
changed the title
|
Ah, my forgotten baby! I got a bit stuck last time trying to find artifacts in test data that would demonstrate usefulness of these reductions. I didn't go very far with this as my main interest was originally in using better estimates of non-gaussian single-trial-power distributions, not necessarily artifact-robustness. But being able to get all the windows without reducing them would be enough gain for me from this PR actually. |
|
Sounds good :) |
|
@mmagnuski you might be interested in this https://github.com/scipy/scipy/pull/8970/files |
6 participants
This PR adds ability to specify reduction to apply to welch windows (#3336)
It's actually part of #3820 with some tests added.
@Eric89GXL I was looking for a branch/stash fixing #3621 and found this instead. But I'll push PR adressing n_fft issue today.
This adds
reductionkeyword argument for reduction across windows, where that can be either'mean','median', callable orNone:'mean'and'median'perform mean or median respectivelyNoneomit the reduction sopsd_welchreturns an array of shapewindows x channels x freqsfor raw andwindows x epochs x channels x freqsfor epochs.scipy.stats.trim_mean)I've added it as
get_reductionto mne.utils. It can also be private if you prefer.I've also added a few tests to
test_psd- they pass locally, but I'll see what CIs have to say.❗️ the windows are the last dimension if a reduction is applied (so if custom function is used it has to reduce across axis=-1), but when reduction is
Nonethey are moved to the -3 dimension. This means that raw output is of diemnsionswindows x channels x freqsand epochs are ofepochs x windows x channels x freqs❗️TODOs:
Closes #3336