對於已經得到的樣本集,核密度估計是一種可以求得樣本的分布的概率密度函數的方法:
通過選取核函數和合適的帶寬,可以得到樣本的distribution probability,在這里核函數選取標准正態分布函數,bandwidth通過AMISE規則選取
具體原理及定義:傳送門 https://en.wikipedia.org/wiki/Density_estimation
MATLAB 代碼實現如下:
1 % Kernel Density Estimation 2 % 只能處理正半軸密度 3 function [t, y_true, tt, y_KDE] = KernelDensityEstimation(x) 4 % clear 5 6 % x = px_last; 7 % x = px_last_tu; 8 %% 9 %參數初始化 10 Max = round(max(x)); %數據中最大值 11 Min = round(min(x)); %數據中最小值 12 Ntotal = length(x); %數據個數 13 tt = 0 : 0.1 : Max; %精確x軸 14 t = 0 : Max; %粗略x軸 15 16 y_KDE = zeros(10 * Max+1, 1); %核密度估計值 17 sum1 = 0; %求和的中間變量 18 %% 19 %計算帶寬h 20 R = 1/(2*sqrt(pi)); 21 m2 = 1; 22 h = 3; 23 % h = (R)^(1/5) / (m2^(2/5) * R^(1/5) * Ntotal^(1/5)); 24 25 %% 26 %計算核密度估計 27 for i = 0 : 0.1 : Max 28 for j = 1 : Ntotal 29 sum1 = sum1 + normpdf(i-x(j)); 30 end 31 y_KDE(round(i*10+1)) = sum1 / (h * Ntotal); 32 sum1 = 0; 33 end 34 35 sum2 = sum(y_KDE)*0.1; %歸一化KDE密度 36 for i = 0 : 0.1 : Max 37 y_KDE(round(i*10+1)) = y_KDE(round(i*10+1))/sum2; 38 end 39 40 %% 41 %計算真實密度的分布 42 y_true = zeros(Max+1,1); 43 for i = 0 : Max 44 for j = 1 : Ntotal 45 if (x(j) < i+1)&&(x(j) >= i) 46 y_true(i+1) = y_true(i+1) + 1; 47 end 48 end 49 y_true(i+1) = y_true(i+1) / Ntotal; 50 end 51 52 %% 53 %繪圖 54 55 % figure(1) %真實密度的分布圖象 56 % bar(t, y_true); 57 % axis([Min Max+1 0 max(y_true)*1.1]); 58 % 59 % figure(2) %核密度估計的密度分布圖象 60 % plot(tt, y_KDE); 61 % axis([Min Max 0 max(y_true)*1.1]);
給定測試數據:
data = [1,2,3,4,5,2,1,2,4,2,1,4,7,4,1,2,4,9,8,7,10,1,2,3,1,0,0,3,6,7,8,9,4]
樣本的條形統計圖和KDE密度分布圖分別如下,可以看到KDE可以較好的還原樣本的分布情況:
真實概率分布圖
KDE密度分布圖