【機率密度函數估計】--- 最大似然估計與Parzen窗函數畫法

【機率密度函數估計】— 最大似然估計與Parsen窗函數畫法

前導知識：【非參數估計—直方圖法、Kn近鄰估計法、Parzen窗法】

1. 最大似然估計

導包

：

import numpy as np
from numpy.linalg import cholesky
import matplotlib.pyplot as plt
import random # 用于随機抽樣
           

設定随機樣本數

：

# 設定随機樣本數
sampleNo = 40;
           

一維資料處理：

mu = np.array([[1]])
Sigma = np.array([[2]])
R = cholesky(Sigma)
s = np.dot(np.random.randn(sampleNo, 1), R) + mu
           

随機抽樣

：

# 随機從40個樣本中抽取20個樣本
n = random.sample(s.tolist(),20)
           

均值估計

：

# 均值估計
u = np.sum(np.array(n))/20
u
           

方差估計

：

# 方差估計
sigma = np.sum((np.array(n)-u)**2)/20
sigma
           

二維資料處理：

mu = np.array([[2, 2]])
Sigma = np.array([[1, 0], [0, 4]])
R = cholesky(Sigma)
s = np.dot(np.random.randn(sampleNo, 2), R) + mu
           

随機抽樣

：

# 随機從40個樣本中抽取20個樣本
n = random.sample(s.tolist(),20)
           

均值向量估計

：

# 均值向量估計
u = np.sum(np.array(n),axis=0)/20
u
           

協方差矩陣估計

：

# 協方差估計
A = np.array(n)-u
B = np.transpose(A)
sigma = np.dot(B,A)/20
sigma
           

array([[0.43709834, 0.32039028],
       [0.32039028, 3.16173429]])
           

2. Parsen窗函數畫法

matlab

實作

主函數：

close all;clear all;
Samples = normrnd(0,1,1,10000); % 從正态分布中産生10000個均值為0，方差為1的樣本
interval = -3:0.01:3; % 劃定橫縱坐标的範圍

index = 1;
for N = [1,10,100]
    for H = [0.25,1,4]
        p = Parsen(Samples,H,N,interval);
        subplot(3,3,index);
        plot(interval,p);
        hold on;
        plot(interval,normpdf(interval,0,1),'r-');
        legend(['h = ',num2str(H),' N = ',num2str(N)]);
        index = index + 1;
    end
end
           

Parsen

函數：

function p = Parsen(Samples,H,N,interval)
% Samples 表示總樣本
% h 表示Parsen視窗大小
% N 是随機采樣的樣本大小（1，10，100）
% x 是密度估計的點

p = zeros(length(interval),1);
h = (H/sqrt(N)); % 半徑

for i = 1 : length(interval)
    b = 0;
    for j = 1 : N
        u = (interval(i) - Samples(j))/h;
        b = b + exp(-u.^2/2)/(sqrt(2*pi)*h); % 一維高斯分布
    end
    p(i) = b / N;
end
end
           

圖示

：

視窗h大小的影響

：

h h h越大，分辨率越低（欠拟合）， h h h越小，穩定性就低些（過拟合）