Index: trunk/examples/Python/lorpe_1d_cv.py =================================================================== --- trunk/examples/Python/lorpe_1d_cv.py (revision 624) +++ trunk/examples/Python/lorpe_1d_cv.py (revision 625) @@ -1,120 +1,120 @@ #!/usr/bin/env python3 """ This example illistrates how to choose LOrPE filter settings in 1-d by cross-validation """ __author__="Igor Volobouev (i.volobouev@ttu.edu)" __version__="1.0" __date__ ="Sep 19 2019" # Perform necessary imports import npstat as ns from npstat_utils import HistoND from npstat_density_estimation_1d import LOrPE1D, FFTKDE1D, \ histoPluginBandwidth1D, BinnedLSCVCalc1D, BinnedRLCVCalc1D import matplotlib.pyplot as plt import numpy as np import math # Some parameters for the script sampleSize = 1000 # Number of points to generate (sample size) xmin = 0.0 # Left edge of the sample histogram xmax = 5.0 # Right edge of the sample histogram numBins = 100 # Number of bins for the sample histogram nDiscrete = 500 # Number of cells for sample discretization. # Should be large enough so that the cell width # is substantially smaller than any physical scale. symbetaPower = -1 # Power parameter of the symmetric beta kernel. # Specify -1 to use a Gaussian kernel instead. pluginDegree = 2 # Degree of polynomial used to generate a simple # plug-in estimate for comparison with cross-validation # Here, we will use the exponential distribution truncated to [xmin, xmax] # in order to generate the random sample d1 = ns.TruncatedDistribution1D(ns.Exponential1D(0.0, 1.0), xmin, xmax) # Generate a random sample according to this density rng = ns.MersenneTwister() sample = d1.generate(rng, sampleSize) # Histogram for displaying the sample h1 = HistoND("int", ns.HistoAxis(numBins, xmin, xmax, "X"), "Sample Histogram", "Counts") any(h1.fill(x, 1) for x in sample) # Histogram for sample discretization h2 = HistoND("int", ns.HistoAxis(nDiscrete, xmin, xmax, "X"), "LOrPE Discretization Histogram", "Counts") any(h2.fill(x, 1) for x in sample) # Object which will perform LOrPE smoother = LOrPE1D(h2) # Add a simple density estimate with a plug-in bandwidth bw0 = histoPluginBandwidth1D(h2, sampleSize, symbetaPower, pluginDegree) print("Plug-in bandwidth is {:.4f} for degree {}".format(bw0, pluginDegree)) yPlugin = smoother.estimateDensity(h2, symbetaPower, pluginDegree, bw0) # Add vanilla KDE for comparison (without mirroring and high order kernels) -kde = FFTKDE1D(h2, True, False) +kde = FFTKDE1D(h2, False, False) bwkde = histoPluginBandwidth1D(h2, sampleSize, symbetaPower, 0) yKDE = kde.estimateDensity(h2, symbetaPower, 0, bwkde) print("KDE bandwidth is {:.4f}".format(bwkde)) # Add a number of filters with different filter degrees and bandwidth # settings for cross-validation calculations for filterDegree in range(1, 7): bw0 = histoPluginBandwidth1D(h2, sampleSize, symbetaPower, filterDegree) print("Building degree {} filters, plug-in bandwidth is {:.4f}".format( filterDegree, bw0)) for bw in np.logspace(math.log10(bw0/4.0), math.log10(bw0*4.0), 50): smoother.addFilter(symbetaPower, filterDegree, bw) # Run the least squares cross-validation (LSCV) lscvCalc = BinnedLSCVCalc1D() bestLSCV, paramsLSCV, yLSCV = smoother.bestFilterByCV(h2, lscvCalc) lscvSymbetaPower, lscvDegree, lscvBandwidth = paramsLSCV print("LSCV: bandwidth is {:.4f}, degree is {}".format(lscvBandwidth, lscvDegree)) # Run the regularized likelihood cross-validation (RLCV) regularisationParameter = 0.5 rlcvCalc = BinnedRLCVCalc1D(regularisationParameter) bestRLCV, paramsRLCV, yRLCV = smoother.bestFilterByCV(h2, rlcvCalc) rlcvSymbetaPower, rlcvDegree, rlcvBandwidth = paramsRLCV print("RLCV: bandwidth is {:.4f}, degree is {}".format(rlcvBandwidth, rlcvDegree)) # Create Matplotlib plot window fig = plt.figure() ax = fig.add_subplot(111) # Plot the sample histogram xh, yh = ns.histoOutline1D(h1) ax.plot(xh, yh, 'k') # Overlay the original density. Scale all densities so that their # plots are compatible with the histogram height. verticalScale = h1.binVolume()*sampleSize xdens, ydens = ns.scanDensity1D(d1, xmin, xmax, nDiscrete+1) ax.plot(xdens, ydens*verticalScale, 'g', label='True density') # Overlay the KDE coords = h2.axis(0).binCenters() ax.plot(coords, yKDE*verticalScale, 'm', label="KDE") # Overlay the LOrPE results with various bandwidth values bandwidths = (bw0, rlcvBandwidth, lscvBandwidth) colors = ('c', 'r', 'b' ) labels = ('Plug-in', 'RLCV', 'LSCV' ) results = (yPlugin, yRLCV, yLSCV ) for bw, color, l, y in zip(bandwidths, colors, labels, results): ax.plot(coords, y*verticalScale, color, label="{} estimate".format(l)) # Add labels and comments plt.xlabel('X') plt.ylabel('Density estimate') plt.title('LOrPE with Cross-Validation') plt.legend() # Display the window plt.show() Index: trunk/examples/Python/kde_1d_by_fft.py =================================================================== --- trunk/examples/Python/kde_1d_by_fft.py (revision 624) +++ trunk/examples/Python/kde_1d_by_fft.py (revision 625) @@ -1,120 +1,120 @@ #!/usr/bin/env python3 """ This example illistrates how to perform constant-bandwidth univariate kernel density estimation (KDE) utilizing fast Fourier transforms (FFT) with the help of the NPStat package. """ __author__="Igor Volobouev (i.volobouev@ttu.edu)" __version__="1.0" __date__ ="Sep 19 2019" # Perform necessary imports import npstat as ns from npstat_utils import HistoND from npstat_density_estimation_1d import FFTKDE1D, histoPluginBandwidth1D, \ BinnedLSCVCalc1D, BinnedRLCVCalc1D import matplotlib.pyplot as plt import numpy as np import math # Some parameters for the script sampleSize = 500 # Number of points to generate (sample size) xmin = -3.5 # Left edge of the sample histogram xmax = 3.5 # Right edge of the sample histogram numBins = 70 # Number of bins for the sample histogram nDiscrete = 1024 # Number of cells for sample discretization. # Should be a power of 2 (this results in the # most efficient application of FFT). symbetaPower = -1 # Power parameter of the symmetric beta kernel. # Specify -1 to use a Gaussian kernel instead. filterDegree = 2 # Degree of polynomial used to generate high-order kernels. # In general, even numbers make the most sense here. # Kernel order will then be filterDegree + 2. The # kernel will be a bona-fide density (i.e., non-negative # everywhere) only if filterDegree is 0. # We will use a 1-d Gassian mixture to generate the random sample. # See the "brute_force_kde_1d.py" script for more comments on this. component1 = ns.GaussianMixtureEntry(3.0, -1.0, 0.5) component2 = ns.GaussianMixtureEntry(7.0, 1.0, 0.5) mixture = ns.GaussianMixture1D(0.0, 1.0, (component1, component2)) # Generate a random sample according to this density rng = ns.MersenneTwister() sample = mixture.generate(rng, sampleSize) # Histogram for displaying the sample h1 = HistoND("int", ns.HistoAxis(numBins, xmin, xmax, "X"), "Sample Histogram", "Counts") any(h1.fill(x, 1) for x in sample) # Histogram for sample discretization. Note that the axis range # is increased, to catch all possible tails. Such an increase may # or may not make sense, depending on application. h2 = HistoND("int", ns.HistoAxis(nDiscrete, 1.5*xmin, 1.5*xmax, "X"), "KDE Discretization Histogram", "Counts") any(h2.fill(x, 1) for x in sample) # Object which will perform KDE smoother = FFTKDE1D(h2) # Make an initial guess for the KDE bandwidth. This guess will use # a robust measure of the sample scale, estimated from the histogram. bw0 = histoPluginBandwidth1D(h2, sampleSize, symbetaPower, filterDegree) -print("Plug-in bandwidth is", bw0) +print("Plug-in bandwidth is {:.4f}".format(bw0)) yPlugin = smoother.estimateDensity(h2, symbetaPower, filterDegree, bw0) # Add a number of kernels with different bandwidth values # for cross-validation calculations for bw in np.logspace(math.log10(bw0/4.0), math.log10(bw0*2.0), 100): smoother.addFilter(symbetaPower, filterDegree, bw) # Run the least squares cross-validation (LSCV) lscvCalc = BinnedLSCVCalc1D() bestLSCV, paramsLSCV, yLSCV = smoother.bestFilterByCV(h2, lscvCalc) lscvSymbetaPower, lscvFilterDegree, lscvBandwidth = paramsLSCV -print("LSCV bandwidth is", lscvBandwidth) +print("LSCV bandwidth is {:.4f}".format(lscvBandwidth)) # Run the regularized likelihood cross-validation (RLCV) regularisationParameter = 0.5 rlcvCalc = BinnedRLCVCalc1D(regularisationParameter) bestRLCV, paramsRLCV, yRLCV = smoother.bestFilterByCV(h2, rlcvCalc) rlcvSymbetaPower, rlcvFilterDegree, rlcvBandwidth = paramsRLCV -print("RLCV bandwidth is", rlcvBandwidth) +print("RLCV bandwidth is {:.4f}".format(rlcvBandwidth)) # We are done with density estimation at this point. The rest of the script # is just plotting results. To begin with, create a Matplotlib plot window. fig = plt.figure() ax = fig.add_subplot(111) # Plot the sample histogram xh, yh = ns.histoOutline1D(h1) ax.plot(xh, yh, 'k') # Overlay the original density. Scale all densities so that their # plots are compatible with the histogram height. verticalScale = h1.binVolume()*sampleSize xdens, ydens = ns.scanDensity1D(mixture, xmin, xmax, nDiscrete+1) ax.plot(xdens, ydens*verticalScale, 'g', label='True density') # Overlay the KDE results with various bandwidth values xkde = h2.axis(0).binCenters() bandwidths = (bw0, rlcvBandwidth, lscvBandwidth) colors = ('c', 'r', 'b' ) labels = ('Plug-in', 'RLCV', 'LSCV' ) results = (yPlugin, yRLCV, yLSCV ) for bw, color, l, y in zip(bandwidths, colors, labels, results): ax.plot(xkde, y*verticalScale, color, label="{} estimate".format(l)) # Set the x range of the plot to that of the sample histogram ax.set_xlim(h1.axis(0).min(), h1.axis(0).max()) # Add labels and comments plt.xlabel('X') plt.ylabel('Density estimate') plt.title('KDE by FFT') plt.legend(loc=2) # Display the window plt.show() Index: trunk/examples/Python/variable_bandwidth_kde_1d.py =================================================================== --- trunk/examples/Python/variable_bandwidth_kde_1d.py (revision 624) +++ trunk/examples/Python/variable_bandwidth_kde_1d.py (revision 625) @@ -1,88 +1,88 @@ #!/usr/bin/env python3 """ This example illistrates how to perform variable-bandwidth univariate kernel density estimation (KDE) with NPStat """ __author__="Igor Volobouev (i.volobouev@ttu.edu)" __version__="1.0" __date__ ="Sep 17 2019" # Perform necessary imports import npstat as ns from npstat_utils import HistoND import matplotlib.pyplot as plt # Some parameters for the script sampleSize = 500 # Number of points to generate (sample size) xmin = -3.5 # Left edge of the sample histogram xmax = 3.5 # Right edge of the sample histogram numBins = 70 # Number of bins for the sample histogram nDiscrete = 1000 # Number of cells for sample discretization symbetaPower = -1 # Power parameter of the symmetric beta kernel. # Specify -1 to use a Gaussian kernel instead. bwFactor = 0.5 # Pilot estimate bandwidth will be set to the product # of the plug-in bandwidth and this factor. The choice # of the factor is empirical here. # We will use a 1-d Gassian mixture to generate the random sample. # See the "brute_force_kde_1d.py" script for more comments on this. component1 = ns.GaussianMixtureEntry(3.0, -1.0, 0.5) component2 = ns.GaussianMixtureEntry(7.0, 1.0, 0.5) mixture = ns.GaussianMixture1D(0.0, 1.0, (component1, component2)) # Generate a random sample according to this density rng = ns.MersenneTwister() sample = mixture.generate(rng, sampleSize) # Histogram for displaying the sample h1 = HistoND("int", ns.HistoAxis(numBins, xmin, xmax, "X"), "Sample Histogram", "Counts") any(h1.fill(x, 1) for x in sample) # Sample discretization for KDE. For smoothing purposes, the bin width # should be substantially smaller than any physical scale of interest. h2 = HistoND("int", ns.HistoAxis(nDiscrete, 2*xmin, 2*xmax, "X"), "KDE Discretization Histogram", "Counts") any(h2.fill(x, 1) for x in sample) # Perform variable-bandwidth KDE. "simpleVariableBandwidthSmooth1D" is the # simplest function to use. It is appropriate when the histogrammed sample # consists of unweighted points. Function "weightedVariableBandwidthSmooth1D" # can be used if the points are weighted (it requires the user to specify # the effective sample size). Function "variableBandwidthSmooth1D" is the # most flexible one but it requires the user to come up with his/her own # pilot bandwidth. Note that Python signatures of these functions are # generated by the "npstat/wrap/variableBandwidthSmooth1DWrap.hh" header # rather than "npstat/stat/variableBandwidthSmooth1D.hh". ykde, bw = ns.simpleVariableBandwidthSmooth1D(h2, symbetaPower, bwFactor) -print("Pilot estimate bandwidth is", bw) +print("Pilot bandwidth is {:.4f}".format(bw)) # Create Matplotlib plot window fig = plt.figure() ax = fig.add_subplot(111) # Plot the sample histogram xh, yh = ns.histoOutline1D(h1) ax.plot(xh, yh, 'k') # Overlay the original density. Scale all densities so that their # plots are compatible with the histogram height. verticalScale = h1.binVolume()*sampleSize xdens, ydens = ns.scanDensity1D(mixture, xmin, xmax, nDiscrete+1) ax.plot(xdens, ydens*verticalScale, 'g', label='True density') # Overlay the KDE result. Scale it to the area of the histogram. ax.plot(h2.axis(0).binCenters(), ykde*verticalScale, 'r', label='KDE') # Set the x range of the plot to that of the sample histogram ax.set_xlim(h1.axis(0).min(), h1.axis(0).max()) # Add labels and comments plt.xlabel('X') plt.ylabel('Density') plt.title('Variable-Bandwidth KDE') plt.legend(loc=2) # Display the window plt.show() Index: trunk/examples/Python/lorpe_1d.py =================================================================== --- trunk/examples/Python/lorpe_1d.py (revision 624) +++ trunk/examples/Python/lorpe_1d.py (revision 625) @@ -1,114 +1,114 @@ #!/usr/bin/env python3 """ This example illistrates how to perform univariate density estimation utilizing local orthogonal polynomial expansion (LOrPE). For details of the method, see the article by D.P. Amali Dassanayake, Igor Volobouev, and A. Alexandre Trindade, "Local orthogonal polynomial expansion for density estimation", Journal of Nonparametric Statistics, 29:4, 806-830 (2017). https://doi.org/10.1080/10485252.2017.1371715 (also arXiv:1505.00275). """ __author__="Igor Volobouev (i.volobouev@ttu.edu)" __version__="1.1" __date__ ="Sep 20 2019" # Perform necessary imports import npstat as ns from npstat_utils import HistoND from npstat_density_estimation_1d import FFTKDE1D, histoPluginBandwidth1D import matplotlib.pyplot as plt # Some parameters for the script sampleSize = 1000 # Number of points to generate (sample size) xmin = 0.0 # Left edge of the sample histogram xmax = 5.0 # Right edge of the sample histogram numBins = 100 # Number of bins for the sample histogram nDiscrete = 1000 # Number of cells for sample discretization symbetaPower = -1 # Power parameter of the symmetric beta weight function # w.r.t. which the orthogonal polynomial system will # be generated. Use -1 to employ a Gaussian instead. filterDegree = 2 # Polynomial degree for the local expansion. This # degree is just a guess here. In general, it can # be determined by cross-validation. bwFactor = 1.0 # LOrPE bandwidth will be set to the product of the # plug-in bandwidth and this factor # Here, we will use the exponential distribution truncated to [xmin, xmax] # in order to generate the random sample d1 = ns.TruncatedDistribution1D(ns.Exponential1D(0.0, 1.0), xmin, xmax) # Generate a random sample according to this density rng = ns.MersenneTwister() sample = d1.generate(rng, sampleSize) # Histogram for displaying the sample h1 = HistoND("int", ns.HistoAxis(numBins, xmin, xmax, "X"), "Sample Histogram", "Counts") any(h1.fill(x, 1) for x in sample) # Sample discretization for LOrPE. For smoothing purposes, the bin width # should be substantially smaller than any physical scale of interest. h2 = HistoND("int", ns.HistoAxis(nDiscrete, xmin, xmax, "X"), "LOrPE Discretization Histogram", "Counts") any(h2.fill(x, 1) for x in sample) # LOrPE boundary handling method. Simple methods (which require no additional # tuning parameters) are: # # "BM_TRUNCATE" -- Just truncate the weight function at the boundary, # as described in the article by Dassanayake et. al. # "BM_STRETCH" -- Stretch the weight function by a factor that smoothly # changes from 1 to 2 as the boundary is approached. # "BM_FOLD" -- Reflect (fold) the weight function at the boundary and # add to the non-reflected part. # "BM_CONSTSQ" -- Stretch the weight function near the boundary in such # a way that the integral of the weight function squared # inside the bounded interval is preserved. # "BM_FOLDSQ" -- Both fold the weight function at the boundary and # preserve the integral of the weight function squared. # "BM_CONSTVAR" -- Stretch the weight function near the boundary so that # the variance of a variable distributed acording to # this weight function remains constant. # "BM_FOLDVAR" -- Both fold the weight function at the boundary and # stretch so that the variance remains constant. # # It is not obvious apriori which boundary handling method will perform # best for any particular sample, so it might be useful to try several. bh = ns.BoundaryHandling("BM_FOLD") # Perform density estimation with LOrPE using a simple high-level # driver function ylorpe, bw = ns.lorpeSmooth1D(h2, symbetaPower, filterDegree, bh, bwFactor) print("LOrPE bandwidth is {:.4f}".format(bw)) # Add vanilla KDE for comparison (without mirroring and high order kernels) -kde = FFTKDE1D(h2, True, False) +kde = FFTKDE1D(h2, False, False) bwkde = bwFactor*histoPluginBandwidth1D(h2, sampleSize, symbetaPower, 0) yKDE = kde.estimateDensity(h2, symbetaPower, 0, bwkde) print("KDE bandwidth is {:.4f}".format(bwkde)) # Plot the sample histogram xh, yh = ns.histoOutline1D(h1) plt.plot(xh, yh, 'k') # Overlay the original density. Scale all densities so that their # plots are compatible with the histogram height. area = h1.integral() xdens, ydens = ns.scanDensity1D(d1, xmin, xmax, nDiscrete+1) plt.plot(xdens, ydens*area, 'g', label='True density') # Overlay the KDE result coords = h2.axis(0).binCenters() plt.plot(coords, yKDE*area, 'b', label='KDE') # Overlay the LOrPE result plt.plot(coords, ylorpe*area, 'r', label='LOrPE') # Add labels and comments plt.xlabel('X') plt.ylabel('Density estimate') plt.title('Density Estimation by KDE and LOrPE') plt.legend() # Display the window plt.show() Index: trunk/examples/Python/brute_force_kde_1d.py =================================================================== --- trunk/examples/Python/brute_force_kde_1d.py (revision 624) +++ trunk/examples/Python/brute_force_kde_1d.py (revision 625) @@ -1,119 +1,119 @@ #!/usr/bin/env python3 """ This example illistrates how to perform constant-bandwidth univariate kernel density estimation (KDE) with a brute-force calculation using the NPStat package. The KDE bandwidth is chosen by cross-validation (CV). This kind of calculation works for relatively small sample sizes as the CV computational complexity is proportional to the sample size squared. """ __author__="Igor Volobouev (i.volobouev@ttu.edu)" __version__="1.1" __date__ ="Sep 15 2019" # Perform necessary imports import npstat as ns from npstat_utils import HistoND from npstat_density_estimation_1d import histoPluginBandwidth1D import matplotlib.pyplot as plt import numpy as np import math # Some parameters for the script sampleSize = 100 # Number of points to generate (sample size) xmin = -3.5 # Left edge of the sample histogram xmax = 3.5 # Right edge of the sample histogram numBins = 35 # Number of bins for the sample histogram scanPoints = 701 # Number of points to use for density scans symbetaPower = 4 # Power parameter of the symmetric beta function used as # the kernel. Use -1 to create a Gaussian kernel instead. filterDegree = 0 # Degree of polynomial used to generate high-order kernels. # In general, even numbers make the most sense here. # Kernel order will then be filterDegree + 2. The # kernel will be a bona-fide density (i.e., non-negative # everywhere) only if filterDegree is 0. # We will use a 1-d Gassian mixture to generate the random sample. The # components of such a mixture are described by the "GaussianMixtureEntry" # class. GaussianMixtureEntry constructor takes three arguments: the # weight of the component, the mean of the Gaussian, and the standard # deviation. Weights can not be negative and the standard deviations # must be positive. The sum of the weights will be normalized internally. component1 = ns.GaussianMixtureEntry(3.0, -1.0, 0.5) component2 = ns.GaussianMixtureEntry(7.0, 1.0, 0.5) mixture = ns.GaussianMixture1D(0.0, 1.0, (component1, component2)) # Generate a random sample according to the Gassian mixture density rng = ns.MersenneTwister() sample = mixture.generate(rng, sampleSize) # Histogram this random sample h1 = HistoND("int", ns.HistoAxis(numBins, xmin, xmax, "X"), "Sample Histogram", "Counts") any(h1.fill(x, 1) for x in sample) # Create the kernel that will be used to perform KDE kernel = ns.KDE1DHOSymbetaKernel(symbetaPower, filterDegree) # Create the object which will be used to estimate the density kde = ns.DoubleKDE1D(kernel, sample) # Make an initial guess for the KDE bandwidth. It will use # a robust measure of the sample scale, estimated from the histogram. bw0 = histoPluginBandwidth1D(h1, sampleSize, symbetaPower, filterDegree) -print("Plug-in bandwidth is", bw0) +print("Plug-in bandwidth is {:.4f}".format(bw0)) # Using the initial bandwidth guess, choose bandwidth by cross-validation. # We can use some numerical optimization method (e.g., scipy.optimize), but # here we will perform just a simple bandwidth scan. Note that for brute-force # KDE the time of cross-validation calculations scales as O(sampleSize^2). bwValues = np.logspace(math.log10(bw0/4.0), math.log10(bw0*2.0), 100) # Run the regularized likelihood cross-validation (RLCV). # For the description of this method, see section 4.2 in # https://doi.org/10.1080/10485252.2017.1371715 regularisationParameter = 0.5 rlcvFunctor = kde.rlcvFunctor(regularisationParameter) rlcv, rlcvBandwidth = max([(rlcvFunctor(bw), bw) for bw in bwValues]) -print("RLCV bandwidth is", rlcvBandwidth) +print("RLCV bandwidth is {:.4f}".format(rlcvBandwidth)) # Run the least squares cross-validation (LSCV). For this type of CV, # the code will need to calculate the integral of density estimate squared. # The parameters of the "lscvFunctor" method specify how this integration # will be performed. The method signature is # # lscvFunctor(xmin, xmax, nSubIntervals, nPoints) # # The interval of integration will be [xmin, xmax]. This interval will # be split into "nSubIntervals" subintervals, and the Gauss-Legendre # quadrature rule with "nPoints" points will be utilized on each subinterval. lscvFunctor = kde.lscvFunctor(2*xmin, 2*xmax, 1, 256) lscv, lscvBandwidth = max([(lscvFunctor(bw), bw) for bw in bwValues]) -print("LSCV bandwidth is", lscvBandwidth) +print("LSCV bandwidth is {:.4f}".format(lscvBandwidth)) # Plot the sample histogram xh, yh = ns.histoOutline1D(h1) plt.plot(xh, yh, 'k') # Overlay the original density. Scale all densities so that their # plots are compatible with the histogram height. verticalScale = h1.binVolume()*sampleSize xdens, ydens = ns.scanDensity1D(mixture, xmin, xmax, scanPoints) plt.plot(xdens, ydens*verticalScale, 'g', label='True density') # Overlay the KDE results with various bandwidth values bandwidths = (bw0, rlcvBandwidth, lscvBandwidth) colors = ('c', 'r', 'b' ) labels = ('Plug-in', 'RLCV', 'LSCV' ) for bw, color, l in zip(bandwidths, colors, labels): x, y = ns.scanKDE1D(kde, xmin, xmax, scanPoints, bw) plt.plot(x, y*verticalScale, color, label="{} estimate".format(l)) # Add labels and comments plt.xlabel('X') plt.ylabel('Density estimate') plt.title('Brute-Force KDE') plt.legend(loc=2) # Display the window plt.show()