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

.. only:: html

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

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

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

.. _sphx_glr_auto_data_analysis_distribution_fitting_plot_smoothing_mixture.py:


Bandwidth sensitivity in kernel smoothing
=========================================

.. GENERATED FROM PYTHON SOURCE LINES 6-41

Introduction
------------

When we have a sample, we may estimate the probability density function of the
underlying distribution using kernel smoothing.
One of the parameters of this method is the bandwidth, which can be either
set by the user, or estimated from the data.
This is especially true when the density is multimodal.
In this example, we consider a bimodal distribution and see how the bandwidth
can change the estimated probability density function.

We consider the distribution:

.. math::
   f_1(x) = w_1 f_A(x) + w_2 f_B(x)


for any :math:`x\in\mathbb{R}` where :math:`f_A` is the density of the Normal distribution
:math:`\mathcal{N}(0,1)`, :math:`f_B` is the density of the Normal distribution
:math:`\mathcal{N}(3/2,(1/3)^2)` and the weights are :math:`w_1 = \frac{3}{4}` and :math:`w_2 = \frac{1}{4}`.

This is a mixture of two Normal distributions: 1/4th of the observations have the :math:`\mathcal{N}(0,1)`
distribution and 3/4th of the observations have the :math:`\mathcal{N}(3/2,(1/3)^2)` distribution.
This example is considered in (Wand, Jones, 1994), page 14, Figure 2.3.

We consider a sample generated from independent realizations of :math:`f_1` and
want to approximate the distribution from kernel smoothing.
More precisely, we want to observe the sensitivity of the resulting density to the
bandwidth.

References
----------

* "Kernel Smoothing", M.P.Wand, M.C.Jones. Chapman and Hall / CRC (1994).


.. GENERATED FROM PYTHON SOURCE LINES 44-53

Generate the mixture by merging two samples
-------------------------------------------

In this section, we show that a mixture of two Normal distributions is
nothing more than the merged sample of two independent Normal distributions.
In order to generate a sample with size :math:`n`, we sample
:math:`\lfloor w_1 n\rfloor` points from the first Normal distribution
:math:`f_A` and :math:`\lfloor w_2 n\rfloor` points from the second Normal
distribution :math:`f_B`. Then we merge the two samples.

.. GENERATED FROM PYTHON SOURCE LINES 55-58

.. code-block:: Python

    import openturns as ot
    import openturns.viewer as otv








.. GENERATED FROM PYTHON SOURCE LINES 59-60

We choose a rather large sample size: :math:`n=1000`.

.. GENERATED FROM PYTHON SOURCE LINES 62-64

.. code-block:: Python

    n = 1000








.. GENERATED FROM PYTHON SOURCE LINES 65-66

Then we define the two Normal distributions and their parameters.

.. GENERATED FROM PYTHON SOURCE LINES 68-73

.. code-block:: Python

    w1 = 0.75
    w2 = 1.0 - w1
    distribution1 = ot.Normal(0.0, 1.0)
    distribution2 = ot.Normal(1.5, 1.0 / 3.0)








.. GENERATED FROM PYTHON SOURCE LINES 74-75

We generate two independent sub-samples from the two Normal distributions.

.. GENERATED FROM PYTHON SOURCE LINES 77-80

.. code-block:: Python

    sample1 = distribution1.getSample(int(w1 * n))
    sample2 = distribution2.getSample(int(w2 * n))








.. GENERATED FROM PYTHON SOURCE LINES 81-82

Then we merge the sub-samples into a larger one with the `add` method of the `Sample` class.

.. GENERATED FROM PYTHON SOURCE LINES 84-88

.. code-block:: Python

    sample = ot.Sample(sample1)
    sample.add(sample2)
    sample.getSize()





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    1000



.. GENERATED FROM PYTHON SOURCE LINES 89-90

In order to see the result, we build a kernel smoothing approximation on the sample. In order to keep it simple, let us use the default bandwidth selection rule.

.. GENERATED FROM PYTHON SOURCE LINES 92-95

.. code-block:: Python

    factory = ot.KernelSmoothing()
    fit = factory.build(sample)








.. GENERATED FROM PYTHON SOURCE LINES 96-99

.. code-block:: Python

    graph = fit.drawPDF()
    view = otv.View(graph)




.. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_smoothing_mixture_001.png
   :alt: plot smoothing mixture
   :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_smoothing_mixture_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 100-103

We see that the distribution of the merged sample has two modes.
However, these modes are not clearly distinct. To distinguish them, we could increase the sample size.
However, it might be interesting to see if the bandwidth selection rule can be better chosen: this is the purpose of the next section.

.. GENERATED FROM PYTHON SOURCE LINES 105-110

Simulation based on a mixture
-----------------------------

Since the distribution that we approximate is a mixture, it will be more convenient to create it from the `Mixture` class.
It takes as input argument a list of distributions and a list of weights.

.. GENERATED FROM PYTHON SOURCE LINES 112-114

.. code-block:: Python

    distribution = ot.Mixture([distribution1, distribution2], [w1, w2])








.. GENERATED FROM PYTHON SOURCE LINES 115-116

Then we generate a sample from it.

.. GENERATED FROM PYTHON SOURCE LINES 118-120

.. code-block:: Python

    sample = distribution.getSample(n)








.. GENERATED FROM PYTHON SOURCE LINES 121-124

.. code-block:: Python

    factory = ot.KernelSmoothing()
    fit = factory.build(sample)








.. GENERATED FROM PYTHON SOURCE LINES 125-127

.. code-block:: Python

    factory.getBandwidth()






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    class=Point name=Unnamed dimension=1 values=[0.208514]
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 128-129

We see that the default bandwidth is close to 0.17.

.. GENERATED FROM PYTHON SOURCE LINES 131-139

.. code-block:: Python

    graph = distribution.drawPDF()
    curve = fit.drawPDF()
    graph.add(curve)
    graph.setColors(["dodgerblue3", "darkorange1"])
    graph.setLegends(["Mixture", "Kernel smoothing"])
    graph.setLegendPosition("upper left")
    view = otv.View(graph)




.. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_smoothing_mixture_002.png
   :alt: plot smoothing mixture
   :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_smoothing_mixture_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 140-141

We see that the result of the kernel smoothing approximation of the mixture is similar to the previous one and could be improved.

.. GENERATED FROM PYTHON SOURCE LINES 143-145

Sensitivity to the bandwidth
----------------------------

.. GENERATED FROM PYTHON SOURCE LINES 147-150

In this section, we observe the sensitivity of the kernel smoothing to the bandwidth.
We consider the three following bandwidths: the small bandwidth 0.05, the large bandwidth 0.54 and 0.18 which is in-between.
For each bandwidth, we use the second optional argument of the `build` method in order to select a specific bandwidth value.

.. GENERATED FROM PYTHON SOURCE LINES 152-176

.. code-block:: Python

    hArray = [0.05, 0.54, 0.18]
    nLen = len(hArray)
    grid = ot.GridLayout(1, len(hArray))
    index = 0
    for i in range(nLen):
        fit = factory.build(sample, [hArray[i]])
        graph = fit.drawPDF()
        exact = distribution.drawPDF()
        curve = exact.getDrawable(0)
        curve.setLegend("Mixture")
        curve.setLineStyle("dashed")
        graph.add(curve)
        graph.setXTitle("X")
        graph.setTitle("h=%.4f" % (hArray[i]))
        graph.setLegends([""])
        graph.setColors(ot.Drawable.BuildDefaultPalette(2))
        bounding_box = graph.getBoundingBox()
        upper_bound = bounding_box.getUpperBound()
        upper_bound[1] = 0.5  # Common y-range
        graph.setBoundingBox(bounding_box)
        grid.setGraph(0, index, graph)
        index += 1

    view = otv.View(grid, figure_kw={"figsize": (10.0, 4.0)})



.. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_smoothing_mixture_003.png
   :alt: , h=0.0500, h=0.5400, h=0.1800
   :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_smoothing_mixture_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 177-180

We see that when the bandwidth is too small, the resulting kernel smoothing has many more modes than the distribution it is supposed to approximate.
When the bandwidth is too large, the approximated distribution is too smooth and has only one mode instead of the expected two modes which are in the mixture distribution.
When the bandwidth is equal to 0.18, the two modes are correctly represented.

.. GENERATED FROM PYTHON SOURCE LINES 182-186

Sensitivity to the bandwidth rule
---------------------------------

The library provides three different rules to compute the bandwidth. In this section, we compare the results that we can get with them.

.. GENERATED FROM PYTHON SOURCE LINES 188-191

.. code-block:: Python

    h1 = factory.computeSilvermanBandwidth(sample)[0]
    h1





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    0.3445636453391276



.. GENERATED FROM PYTHON SOURCE LINES 192-195

.. code-block:: Python

    h2 = factory.computePluginBandwidth(sample)[0]
    h2





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    0.2021709523195656



.. GENERATED FROM PYTHON SOURCE LINES 196-199

.. code-block:: Python

    h3 = factory.computeMixedBandwidth(sample)[0]
    h3





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    0.20851397168332242



.. GENERATED FROM PYTHON SOURCE LINES 200-202

.. code-block:: Python

    factory.getBandwidth()[0]





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    0.18



.. GENERATED FROM PYTHON SOURCE LINES 203-207

We see that the default rule is the "Mixed" rule. This is because the sample is in dimension 1 and the sample size is quite large.
For a small sample in 1 dimension, the "Plugin" rule would have been used.

The following script compares the results produced by the three rules.

.. GENERATED FROM PYTHON SOURCE LINES 209-238

.. code-block:: Python

    hArray = [h1, h2, h3]
    legends = ["Silverman", "Plugin", "Mixed"]
    nLen = len(hArray)
    grid = ot.GridLayout(1, len(hArray))
    index = 0
    for i in range(nLen):
        fit = factory.build(sample, [hArray[i]])
        graph = fit.drawPDF()
        exact = distribution.drawPDF()
        curve = exact.getDrawable(0)
        curve.setLegend("Mixture")
        curve.setLineStyle("dashed")
        graph.add(curve)
        graph.setLegends([""])
        graph.setTitle("h=%.4f, %s" % (hArray[i], legends[i]))
        graph.setXTitle("X")
        graph.setColors(ot.Drawable.BuildDefaultPalette(2))
        if i > 0:
            graph.setYTitle("")
        bounding_box = graph.getBoundingBox()
        upper_bound = bounding_box.getUpperBound()
        upper_bound[1] = 0.5  # Common y-range
        graph.setBoundingBox(bounding_box)
        grid.setGraph(0, index, graph)
        index += 1

    view = otv.View(grid, figure_kw={"figsize": (10.0, 4.0)})

    otv.View.ShowAll()



.. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_smoothing_mixture_004.png
   :alt: , h=0.3446, Silverman, h=0.2022, Plugin, h=0.2085, Mixed
   :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_smoothing_mixture_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 239-241

We see that the bandwidth produced by Silverman's rule is too large, leading to an oversmoothed distribution.
The results produced by the Plugin and Mixed rules are comparable in this case.


.. _sphx_glr_download_auto_data_analysis_distribution_fitting_plot_smoothing_mixture.py:

.. only:: html

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

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

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

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

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