.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples/plot_epsilon_selection.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_examples_plot_epsilon_selection.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_plot_epsilon_selection.py:


Choosing epsilon with cross-validation
--------------------------------------

The ``epsilon`` parameter of :meth:`pysecs.SECS.fit` controls how many
singular values are kept in the inversion: too small and the fit chases
noise, too large and real structure is smoothed away. Leave-one-out
cross-validation gives an objective way to choose it: hold out one
station, fit the rest, and measure the prediction error at the held-out
station.

.. GENERATED FROM PYTHON SOURCE LINES 12-63



.. image-sg:: /examples/images/sphx_glr_plot_epsilon_selection_001.png
   :alt: Cross-validated choice of the regularization parameter
   :srcset: /examples/images/sphx_glr_plot_epsilon_selection_001.png
   :class: sphx-glr-single-img





.. code-block:: Python


    import matplotlib.pyplot as plt
    import numpy as np

    from pysecs import SECS, T_df


    R_E = 6378e3
    rng = np.random.default_rng(42)

    # A grid of divergence-free SECs at 110 km altitude
    lats = np.linspace(35, 55, 8)
    lons = np.linspace(-10, 20, 10)
    lat_g, lon_g = np.meshgrid(lats, lons, indexing="ij")
    sec_locs = np.column_stack(
        [lat_g.ravel(), lon_g.ravel(), np.full(lat_g.size, R_E + 110e3)]
    )

    # Synthetic "true" current system and noisy ground observations
    true_amps = rng.normal(0, 1e4, sec_locs.shape[0])
    nobs = 25
    obs_locs = np.column_stack(
        [rng.uniform(36, 54, nobs), rng.uniform(-9, 19, nobs), np.full(nobs, R_E)]
    )
    B_true = np.tensordot(true_amps, T_df(obs_locs, sec_locs), (0, 2))
    noise = 0.05 * np.max(np.abs(B_true))
    B_obs = B_true + rng.normal(0, noise, B_true.shape)

    # Leave-one-out cross-validation over a range of epsilons
    epsilons = np.logspace(-4, 0, 20)
    cv_error = np.zeros_like(epsilons)
    for i, epsilon in enumerate(epsilons):
        for k in range(nobs):
            keep = np.arange(nobs) != k
            secs = SECS(sec_df_loc=sec_locs)
            secs.fit(obs_locs[keep], B_obs[keep], epsilon=epsilon)
            pred = secs.predict_B(obs_locs[[k]])
            cv_error[i] += np.sum((pred - B_obs[k]) ** 2)
    cv_error = np.sqrt(cv_error / (nobs * 3))

    best = epsilons[np.argmin(cv_error)]

    fig, ax = plt.subplots()
    ax.loglog(epsilons, cv_error * 1e9, marker="o")
    ax.axvline(best, color="tab:red", ls="--", label=f"best epsilon = {best:.2g}")
    ax.axhline(noise * 1e9, color="gray", ls=":", label="noise level")
    ax.set_xlabel("epsilon")
    ax.set_ylabel("leave-one-out RMSE (nT)")
    ax.set_title("Cross-validated choice of the regularization parameter")
    ax.legend()
    plt.show()


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 1.060 seconds)


.. _sphx_glr_download_examples_plot_epsilon_selection.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_epsilon_selection.ipynb <plot_epsilon_selection.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_epsilon_selection.py <plot_epsilon_selection.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_epsilon_selection.zip <plot_epsilon_selection.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_