Non-parametric density estimation for a bimodal sample

kernel_density_estimation.py

# -*- coding: utf-8 -*-
"""
Created on Tue Nov 17 20:58:21 2020

@author: Localuser
"""


import matplotlib.pyplot as plt
import numpy as np
from sklearn.neighbors import KernelDensity
import math

#creating and plotting a bimodal data sample
sample1 = np.random.normal(loc = 20, scale = 5, size = 300)
sample2 = np.random.normal(loc = 40, scale = 5, size = 700)
sample = np.hstack((sample1, sample2))
plt.hist(sample, bins = 50)
plt.show()

# fitting the density
model = KernelDensity(bandwidth = 2, kernel = 'gaussian')
sample = sample.reshape((len(sample), 1))
model.fit(sample)

# sample probabilities for a range of outcomes
values = np.asarray([value for value in range(1, 60)])
values = values.reshape((len(values), 1))
probabilities = model.score_samples(values)
probabilities = math.exp(probabilities)

# plot the histogram and pdf
plt.hist(sample, bins=50, density=True)
plt.plot(values[:], probabilities)
plt.show()