天天看點

高斯估計python和matlab

高斯估計 我也忘了具體是啥了

import numpy as np
import math
NDATA=3000
MU=[1,2]
COV=[[0.25,0],[0,0.64]]
SIGMA=[0.5,0.8]
#data=np.random.multivariate_normal(MU,SIGMA,NDATA)
#print(data)
E=np.zeros((NDATA,1,2))
DATA=np.zeros((NDATA,1,2))
COV_R=np.zeros((NDATA,2,2))
MU_R=0
#print(DATA)
for i in range(NDATA):
  DATA[i] = np.random.multivariate_normal(MU, COV)
for i in range(NDATA):
 MU_R=MU_R+DATA[i,0]
MU_R=MU_R/NDATA    #抽取值
for i in range(NDATA):
  #print(DATA[i])
  E[i]=DATA[i]-MU_R #E[i]為第i個樣本的偏差 [X-MU,Y-MU]
  #print(E[i])
for i in range(NDATA):
 #print(E[i])
 COV_R[i]=E[i]*E[i].T
 #print(COV_R[i])
COV_1=0
COV_2=0
for i in range(NDATA):  #均值用mean() 然我沒看明白
 COV_1=COV_R[i][0][0]+COV_1
 #print(COV_R[i][0][0])
 COV_2=COV_R[i][1][1]+COV_2
COV_1=math.sqrt(COV_1/NDATA)
COV_2=math.sqrt(COV_2/NDATA)
print(MU_R)
print(COV_1)
print(COV_2)
           

下面是老師給的matlab版,我相當于複寫吧,拿python寫了一遍。

% 本示例用于最大似然估計的教學示例
% 描述:1)生成過程:根據給定參數的正态分布生成二維資料點樣本
% ···2)估計過程:根據最大似然估計算法估計參數值
% 作者:王家寶

% 設定随機數種子,使程式可重複 10次産生的随機數相同
rng(10);

% 生成過程
NDATA = 10; % 資料點數目
MU = [1,2]'; % 分布均值 '代表轉置
COV = [0.5,0;0,0.8]; % 協方差:假設各次元互相獨立
SIGMA = diag(COV); % 方差:假設各次元互相獨立 提取主對角線
DATA = zeros(2,NDATA); %DATA為2*NDATA的零矩陣
for k = 1 : NDATA
    % 依次生成資料
    DATA(:,k) = normrnd(MU,SIGMA);%生成服從正态分布的随機數
end
disp('均值:');
disp(MU);
% disp('協方差:');
% disp(COV);
disp('方差:');
disp(SIGMA);

scatter(DATA(1,:),DATA(2,:)); % 繪制樣本點 第一行為橫軸 第二行為縱軸
box on; %右邊上邊有邊框了
axis equal;%坐标軸機關相同

% 估計過程
MU_EST = mean(DATA,2); % 計算均值 計算每行均值 MU_EST為估計均值 DATA為資料
DATA_DIF = bsxfun(@minus,DATA,MU_EST); %相減為資料與均值的內插補點
COV_SQUARE = zeros(2,2,NDATA);%2*2*NDATDA的矩陣
for k = 1 : NDATA
    % 計算每個樣本的協方差
    COV_SQUARE(:,:,k) = DATA_DIF(:,k)*DATA_DIF(:,k)';
end
COV_SQUARE_EST = mean(COV_SQUARE,3);%對NDATA個樣本取均值 c11:E(xi-MUx)^2 c22:E(yi-MUy)^2
SIGMA_SQUARE_EST = diag(COV_SQUARE_EST);%取對角線
disp('估計均值:');
disp(MU_EST);
% disp('估計協方差:');
% disp(sqrt(COV_SQUARE_EST));
disp('估計方差:');
disp(sqrt(SIGMA_SQUARE_EST));

% 注:系統自帶的參數估計,假設資料各次元互相獨立
[mu,s] = normfit(DATA');