.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/noise_alpha/plot_simulate_alphacsc.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_noise_alpha_plot_simulate_alphacsc.py: =========================== Alpha CSC on simulated data =========================== This example demonstrates alphaCSC [1]_ on simulated data. .. [1] Jas, M., Dupré La Tour, T., Şimşekli, U., & Gramfort, A. (2017). `Learning the Morphology of Brain Signals Using Alpha-Stable Convolutional Sparse Coding `_. Advances in Neural Information Processing Systems (NIPS), pages 1099--1108. .. GENERATED FROM PYTHON SOURCE LINES 14-22 .. code-block:: Python # Authors: Mainak Jas # Tom Dupre La Tour # Umut Simsekli # Alexandre Gramfort # # License: BSD (3-clause) .. GENERATED FROM PYTHON SOURCE LINES 23-24 Let us first define the parameters of our model. .. GENERATED FROM PYTHON SOURCE LINES 24-33 .. code-block:: Python n_times_atom = 64 # L n_times = 512 # T n_atoms = 2 # K n_trials = 100 # N alpha = 1.2 reg = 0.1 .. GENERATED FROM PYTHON SOURCE LINES 34-35 Next, we define the parameters for alpha CSC .. GENERATED FROM PYTHON SOURCE LINES 35-41 .. code-block:: Python n_iter_global = 10 n_iter_optim = 50 n_iter_mcmc = 100 n_burnin_mcmc = 50 .. GENERATED FROM PYTHON SOURCE LINES 42-43 Here, we simulate the data .. GENERATED FROM PYTHON SOURCE LINES 43-50 .. code-block:: Python from alphacsc.simulate import simulate_data random_state_simulate = 1 X, ds_true, z_true = simulate_data(n_trials, n_times, n_times_atom, n_atoms, random_state_simulate) .. GENERATED FROM PYTHON SOURCE LINES 51-52 Add some noise and corrupt some trials even with impulsive noise .. GENERATED FROM PYTHON SOURCE LINES 52-71 .. code-block:: Python from scipy.stats import levy_stable from alphacsc import check_random_state fraction_corrupted = 0.02 n_corrupted_trials = int(fraction_corrupted * n_trials) # Add stationary noise: rng = check_random_state(random_state_simulate) X += 0.01 * rng.randn(*X.shape) noise_level = 0.005 # add impulsive noise idx_corrupted = rng.randint(0, n_trials, size=n_corrupted_trials) X[idx_corrupted] += levy_stable.rvs(alpha, 0, loc=0, scale=noise_level, size=(n_corrupted_trials, n_times), random_state=random_state_simulate) .. GENERATED FROM PYTHON SOURCE LINES 72-73 and then alpha CSC on the same data .. GENERATED FROM PYTHON SOURCE LINES 73-83 .. code-block:: Python from alphacsc import learn_d_z_weighted d_hat, z_hat, Tau = learn_d_z_weighted( X, n_atoms, n_times_atom, reg=reg, alpha=alpha, solver_d_kwargs=dict(factr=100), n_iter_global=n_iter_global, n_iter_optim=n_iter_optim, init_tau=True, n_iter_mcmc=n_iter_mcmc, n_burnin_mcmc=n_burnin_mcmc, random_state=60, n_jobs=1, verbose=1) .. rst-class:: sphx-glr-script-out .. code-block:: none Global Iter: 0/10 V_0/50 /github/workspace/alphacsc/update_d.py:172: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.) lambd_hats[k] = lambd_hat ................................................. Global Iter: 1/10 V_0/50 ................................................. Global Iter: 2/10 V_0/50 ................................................. Global Iter: 3/10 V_0/50 ................................................. Global Iter: 4/10 V_0/50 ................................................. Global Iter: 5/10 V_0/50 ................................................. Global Iter: 6/10 V_0/50 ................................................. Global Iter: 7/10 V_0/50 ................................................. Global Iter: 8/10 V_0/50 ................................................. Global Iter: 9/10 V_0/50 ................................................. .. GENERATED FROM PYTHON SOURCE LINES 84-86 Finally, let's compare the results. Now, it works even in the presence of impulsive noise. .. GENERATED FROM PYTHON SOURCE LINES 86-96 .. code-block:: Python import matplotlib.pyplot as plt plt.figure() plt.plot(d_hat.T, 'b', label=r'$\alpha$CSC') plt.plot(ds_true.T, 'k--', label='True atoms') handles, labels = plt.gca().get_legend_handles_labels() plt.legend(handles[::2], labels[::2], loc='best') plt.show() .. image-sg:: /auto_examples/noise_alpha/images/sphx_glr_plot_simulate_alphacsc_001.png :alt: plot simulate alphacsc :srcset: /auto_examples/noise_alpha/images/sphx_glr_plot_simulate_alphacsc_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 97-99 We can even visualize the weights to see what time points were downweighted by the algorithm .. GENERATED FROM PYTHON SOURCE LINES 99-113 .. code-block:: Python fig, axes = plt.subplots(2, 1, sharex=True) axes[0].set_xlim([-20, n_times]) axes[0].set_ylim([0, 2]) axes[1].set_ylim([-0.5, 0.3]) for t in [2, 250]: axes[0].axvline(t, linestyle='-.') axes[1].axvline(t, linestyle='-.') axes[0].plot(1 / Tau[idx_corrupted[0], :]) axes[0].set_ylabel('1 / weights') axes[1].plot(X[idx_corrupted[0], :]) axes[1].set_ylabel('Corrupted trial') .. image-sg:: /auto_examples/noise_alpha/images/sphx_glr_plot_simulate_alphacsc_002.png :alt: plot simulate alphacsc :srcset: /auto_examples/noise_alpha/images/sphx_glr_plot_simulate_alphacsc_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(31.222222222222214, 0.5, 'Corrupted trial') .. rst-class:: sphx-glr-timing **Total running time of the script:** (3 minutes 27.496 seconds) .. _sphx_glr_download_auto_examples_noise_alpha_plot_simulate_alphacsc.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_simulate_alphacsc.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_simulate_alphacsc.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_simulate_alphacsc.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_