How to compute the probability (e.g. 5%, 10%, 90%) of the Kernel density function?

user_19523905 0
Viewed 106

I attempted to plot the kernel density distribution (Gaussian) curve along with the histogram plot of two data set in python.

However, in my script the estimation of 95% (data1: marked by red color vertical line) and 5% (data2: marked by black color vertical line) is very time-consuming, e.g. I need to test different limits [detail explanation in code, where I need to change the upper limited] to get the 95% and 5% probability of the kernel density curve.

May someone help out me here and suggest possible way out fixed this issue or another approach to plot the kernel density curve along with its 95% and 5% probability.

Thank you!

My script is here.

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from sklearn.neighbors import KernelDensity
from scipy import stats
data1 = result['95_24'] # data 1
data2 = result['5_24'] # data 2 
def plot_prob_density(data1, data2, x_start1, x_end1):
    fig, (ax1) = plt.subplots(1, 1, figsize=(6,5), sharey=False)

    unit = 1.5
    x = np.linspace(-20, 20, 1000)[:, np.newaxis]
    
    # Hisogram plot of data 
    ax1.hist(data1, bins=np.linspace(-20,20,40), density=True, color='r', alpha=0.4)
    ax1.hist(data2, bins=np.linspace(-20,20,40), density=True, color='k', alpha=0.4)


    # kernel density estimation
    kd_data1 = KernelDensity(kernel='gaussian', bandwidth=1.8).fit(data1)
    kd_data2 = KernelDensity(kernel='gaussian', bandwidth=1.8).fit(data2)

    kd_vals_data1 = np.exp(kd_data1.score_samples(x))
    kd_vals_data2 = np.exp(kd_data2.score_samples(x))
    
    # density plot 
    ax1.plot(x, kd_vals_data1, color='r', label='$Na$', linewidth=2)
    ax1.plot(x, kd_vals_data2, color='k', label='$Λ$', linewidth = 2) 
  
    #  using the function get probability)
    ax1.axvline(x=x_end1,color='red',linestyle='dashed', linewidth = 3, label='$β_{95\%}$')
    ax1.axvline(x=x_start1,color='k',linestyle='dashed', linewidth = 3, label='$β_{5\%}$')
    

    # Show the plots
    ax1.set_ylabel('Probability density', fontsize=12)
    ax1.set_xlabel('Beta', fontsize=12)
    ax1.set_xlim([-20, 20])
    ax1.set_ylim(0, 0.3)
    ax1.set_yticks([0, 0.1, 0.2, 0.3]) 
    ax1.set_xticks([-20, 20, -10, 10, 0])
    
    ax1.legend(fontsize=12, loc='upper left', frameon=False)
    
    fig.tight_layout()
    gc.collect()
    return kd_data1, kd_data2,

# Calculation of 95% and 5 % for data1 and data2 Kernel density curve


def get_probability(start_value, end_value, eval_points, kd):
    
    # Number of evaluation points 
    N = eval_points                                      
    step = (end_value - start_value) / (N - 1)  # Step size

    x = np.linspace(start_value, end_value, N)[:, np.newaxis]  # Generate values in the range
    kd_vals = np.exp(kd.score_samples(x))  # Get PDF values for each x
    probability = np.sum(kd_vals * step)  # Approximate the integral of the PDF
    return probability.round(4)

data1 = np.array(data1).reshape(-1, 1)
data2 = np.array(data2).reshape(-1, 1)


kd_data1, kd_data2= plot_prob_density(data1, data2, x_start1=-2.2, x_end1=5.3)



# ##############################
print('Beta-95%: {}'
      .format(get_probability(start_value = -20, 
                              end_value = 5.3, 
                              eval_points = 1000, 
                              kd = kd_data1)))
# here, I modify the end-value every time and then see teh output  #value, when it reached to 95% then i took taht values as 95% #confidence, however this is very confsuing, i want to compute this 95% directly and same for 5% probbaility, computed below:
print('Beta-5%: {}\n'
      .format(get_probability(start_value = -20, 
                              end_value = -2.2, 
                              eval_points = 1000, 
                              kd = kd_data2)))
####################################################################

plt.savefig("Ev_test.png")

The pictorial representation is also attached here. Histogram and kernel density plot along with its 95% and 5% probability limits highlighted with red and black vertical bold lines:

Figure 1

pandas probability python
1 Answers

Here is the possible way out to fix this issue. Additionally, the proposed method it has error in percentile calculation, therefore i recommend not to use that:

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import gaussian_kde
import seaborn as sns
from sklearn.neighbors import KernelDensity
%matplotlib inline
import numpy as np
from scipy import stats
import statsmodels.api as sm
import matplotlib.pyplot as plt
from statsmodels.distributions.mixture_rvs import mixture_rvs
from scipy.stats import norm
import numpy as np



fig = plt.figure(figsize=(4, 4), dpi=300)
ax = fig.add_subplot(111)

# Plot the histogram
ax.hist(data8,bins=20,zorder=1,color="r",density=True,alpha=0.6,)
ax.hist(data7,bins=20,zorder=1,color="black",density=True,alpha=0.6,)

# kde.fit()
kde = sm.nonparametric.KDEUnivariate(data8)
kde1 = sm.nonparametric.KDEUnivariate(data7)

# Plot the KDE for various bandwidths
for bandwidth in [1.8]:
    kde.fit(bw=bandwidth)
    kde1.fit(bw=bandwidth)# Estimate the densities
    ax.plot(kde.support, kde.density,"-",lw=2,color="r",zorder=10, alpha=0.6, label="Data1")
    ax.plot(kde1.support, kde1.density,"-",lw=2,color="black",zorder=10, alpha=0.6, label="Data2")


ax.legend(loc="best")
ax.set_xlim([-20, 40])
ax.set_ylim([0, 0.3])
ax.grid(False)

# Probabilities calculation 

quantiles_mesh = np.linspace(0,1,len(kde.density))
fig = plt.figure(figsize=(2, 2), dpi=300)
plt.plot(quantiles_mesh, kde.icdf)

data_1_95=  np.percentile(kde1.icdf, 95)
data_2_5=  np.percentile(kde2.icdf, 5)
ax.axvline(x=data_1_95,color='red',linestyle='dashed', linewidth = 2)
ax.axvline(x=data_2_5,color='k',linestyle='dashed', linewidth = 2)
#plt.savefig("KDE_Plot.png")
user_195239050
Related