{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# Bandwidth sensitivity in kernel smoothing\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction\n\nWhen we have a sample, we may estimate the probability density function of the\nunderlying distribution using kernel smoothing.\nOne of the parameters of this method is the bandwidth, which can be either\nset by the user, or estimated from the data.\nThis is especially true when the density is multimodal.\nIn this example, we consider a bimodal distribution and see how the bandwidth\ncan change the estimated probability density function.\n\nWe consider the distribution:\n\n\\begin{align}f_1(x) = w_1 f_A(x) + w_2 f_B(x)\\end{align}\n\n\nfor any $x\\in\\mathbb{R}$ where $f_A$ is the density of the Normal distribution\n$\\mathcal{N}(0,1)$, $f_B$ is the density of the Normal distribution\n$\\mathcal{N}(3/2,(1/3)^2)$ and the weights are $w_1 = \\frac{3}{4}$ and $w_2 = \\frac{1}{4}$.\n\nThis is a mixture of two Normal distributions: 1/4th of the observations have the $\\mathcal{N}(0,1)$\ndistribution and 3/4th of the observations have the $\\mathcal{N}(3/2,(1/3)^2)$ distribution.\nThis example is considered in (Wand, Jones, 1994), page 14, Figure 2.3.\n\nWe consider a sample generated from independent realizations of $f_1$ and\nwant to approximate the distribution from kernel smoothing.\nMore precisely, we want to observe the sensitivity of the resulting density to the\nbandwidth.\n\n## References\n\n* \"Kernel Smoothing\", M.P.Wand, M.C.Jones. Chapman and Hall / CRC (1994).\n\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Generate the mixture by merging two samples\n\nIn this section, we show that a mixture of two Normal distributions is\nnothing more than the merged sample of two independent Normal distributions.\nIn order to generate a sample with size $n$, we sample\n$\\lfloor w_1 n\\rfloor$ points from the first Normal distribution\n$f_A$ and $\\lfloor w_2 n\\rfloor$ points from the second Normal\ndistribution $f_B$. Then we merge the two samples.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import openturns as ot\nimport openturns.viewer as otv"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We choose a rather large sample size: $n=1000$.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"n = 1000"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then we define the two Normal distributions and their parameters.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"w1 = 0.75\nw2 = 1.0 - w1\ndistribution1 = ot.Normal(0.0, 1.0)\ndistribution2 = ot.Normal(1.5, 1.0 / 3.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We generate two independent sub-samples from the two Normal distributions.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"sample1 = distribution1.getSample(int(w1 * n))\nsample2 = distribution2.getSample(int(w2 * n))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then we merge the sub-samples into a larger one with the `add` method of the `Sample` class.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"sample = ot.Sample(sample1)\nsample.add(sample2)\nsample.getSize()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"factory = ot.KernelSmoothing()\nfit = factory.build(sample)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"graph = fit.drawPDF()\nview = otv.View(graph)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that the distribution of the merged sample has two modes.\nHowever, these modes are not clearly distinct. To distinguish them, we could increase the sample size.\nHowever, it might be interesting to see if the bandwidth selection rule can be better chosen: this is the purpose of the next section.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simulation based on a mixture\n\nSince the distribution that we approximate is a mixture, it will be more convenient to create it from the `Mixture` class.\nIt takes as input argument a list of distributions and a list of weights.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"distribution = ot.Mixture([distribution1, distribution2], [w1, w2])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then we generate a sample from it.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"sample = distribution.getSample(n)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"factory = ot.KernelSmoothing()\nfit = factory.build(sample)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"factory.getBandwidth()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that the default bandwidth is close to 0.17.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"graph = distribution.drawPDF()\ncurve = fit.drawPDF()\ngraph.add(curve)\ngraph.setColors([\"dodgerblue3\", \"darkorange1\"])\ngraph.setLegends([\"Mixture\", \"Kernel smoothing\"])\ngraph.setLegendPosition(\"topleft\")\nview = otv.View(graph)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that the result of the kernel smoothing approximation of the mixture is similar to the previous one and could be improved.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Sensitivity to the bandwidth\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this section, we observe the sensitivity of the kernel smoothing to the bandwidth.\nWe consider the three following bandwidths: the small bandwidth 0.05, the large bandwidth 0.54 and 0.18 which is in-between.\nFor each bandwidth, we use the second optional argument of the `build` method in order to select a specific bandwidth value.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"hArray = [0.05, 0.54, 0.18]\nnLen = len(hArray)\ngrid = ot.GridLayout(1, len(hArray))\nindex = 0\nfor i in range(nLen):\n fit = factory.build(sample, [hArray[i]])\n graph = fit.drawPDF()\n exact = distribution.drawPDF()\n curve = exact.getDrawable(0)\n curve.setLegend(\"Mixture\")\n curve.setLineStyle(\"dashed\")\n graph.add(curve)\n graph.setXTitle(\"X\")\n graph.setTitle(\"h=%.4f\" % (hArray[i]))\n graph.setLegends([\"\"])\n graph.setColors(ot.Drawable.BuildDefaultPalette(2))\n bounding_box = graph.getBoundingBox()\n upper_bound = bounding_box.getUpperBound()\n upper_bound[1] = 0.5 # Common y-range\n graph.setBoundingBox(bounding_box)\n grid.setGraph(0, index, graph)\n index += 1\n\nview = otv.View(grid, figure_kw={\"figsize\": (10.0, 4.0)})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\nWhen 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.\nWhen the bandwidth is equal to 0.18, the two modes are correctly represented.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Sensitivity to the bandwidth rule\n\nThe library provides three different rules to compute the bandwidth. In this section, we compare the results that we can get with them.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"h1 = factory.computeSilvermanBandwidth(sample)[0]\nh1"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"h2 = factory.computePluginBandwidth(sample)[0]\nh2"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"h3 = factory.computeMixedBandwidth(sample)[0]\nh3"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"factory.getBandwidth()[0]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\nFor a small sample in 1 dimension, the \"Plugin\" rule would have been used.\n\nThe following script compares the results produced by the three rules.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"hArray = [h1, h2, h3]\nlegends = [\"Silverman\", \"Plugin\", \"Mixed\"]\nnLen = len(hArray)\ngrid = ot.GridLayout(1, len(hArray))\nindex = 0\nfor i in range(nLen):\n fit = factory.build(sample, [hArray[i]])\n graph = fit.drawPDF()\n exact = distribution.drawPDF()\n curve = exact.getDrawable(0)\n curve.setLegend(\"Mixture\")\n curve.setLineStyle(\"dashed\")\n graph.add(curve)\n graph.setLegends([\"\"])\n graph.setTitle(\"h=%.4f, %s\" % (hArray[i], legends[i]))\n graph.setXTitle(\"X\")\n graph.setColors(ot.Drawable.BuildDefaultPalette(2))\n if i > 0:\n graph.setYTitle(\"\")\n bounding_box = graph.getBoundingBox()\n upper_bound = bounding_box.getUpperBound()\n upper_bound[1] = 0.5 # Common y-range\n graph.setBoundingBox(bounding_box)\n grid.setGraph(0, index, graph)\n index += 1\n\nview = otv.View(grid, figure_kw={\"figsize\": (10.0, 4.0)})\n\notv.View.ShowAll()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that the bandwidth produced by Silverman's rule is too large, leading to an oversmoothed distribution.\nThe results produced by the Plugin and Mixed rules are comparable in this case.\n\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.6"
}
},
"nbformat": 4,
"nbformat_minor": 0
}