from functools import partial import numpy as np from scipy import stats import matplotlib.pyplot as plt def my_kde_bandwidth(obj, fac=1./5): """We use Scott's Rule, multiplied by a constant factor.""" return np.power(obj.n, -1./(obj.d+4)) * fac loc1, scale1, size1 = (-2, 1, 175) loc2, scale2, size2 = (2, 0.2, 50) x2 = np.concatenate([np.random.normal(loc=loc1, scale=scale1, size=size1), np.random.normal(loc=loc2, scale=scale2, size=size2)]) x_eval = np.linspace(x2.min() - 1, x2.max() + 1, 500) kde = stats.gaussian_kde(x2) kde2 = stats.gaussian_kde(x2, bw_method='silverman') kde3 = stats.gaussian_kde(x2, bw_method=partial(my_kde_bandwidth, fac=0.2)) kde4 = stats.gaussian_kde(x2, bw_method=partial(my_kde_bandwidth, fac=0.5)) pdf = stats.norm.pdf bimodal_pdf = pdf(x_eval, loc=loc1, scale=scale1) * float(size1) / x2.size + \ pdf(x_eval, loc=loc2, scale=scale2) * float(size2) / x2.size fig = plt.figure(figsize=(8, 6)) ax = fig.add_subplot(111) ax.plot(x2, np.zeros(x2.shape), 'b+', ms=12) ax.plot(x_eval, kde(x_eval), 'k-', label="Scott's Rule") ax.plot(x_eval, kde2(x_eval), 'b-', label="Silverman's Rule") ax.plot(x_eval, kde3(x_eval), 'g-', label="Scott * 0.2") ax.plot(x_eval, kde4(x_eval), 'c-', label="Scott * 0.5") ax.plot(x_eval, bimodal_pdf, 'r--', label="Actual PDF") ax.set_xlim([x_eval.min(), x_eval.max()]) ax.legend(loc=2) ax.set_xlabel('x') ax.set_ylabel('Density') plt.show()