计算数据的经验分布函数与MATLAB作图

计算数据的经验分布函数与MATLAB作图

写在前面

因为某些原因，需要处理某些数据，比如说某项测量数据与理论值的误差，我们就需要检验误差是否符合正态分布。最直观的方法就是直接做出经验分布函数的图来进行观察。所以这里简单的写一下如何做出经验分布图。

第一步，做经验分布图

对数据按升序排列

首先我们将所有一共n个数据按照升序排列。之所以升序，是因为这样方便我们后续计算。

将区间[a, b]分成小段

a是所有数据的最小值，b是所有数据的最大值，我们将所有数据分成小区间，比如说将[a,b]分为

1. [a,a+(b-a)/100]
2. [a, a+2*(b-a)/100]
3. …
4. [a, b]

注意，区间左值永远是a，然后计算每个区间内落入数据点的个数，这样就得到了一个类似于频数分布直方图的东西。我们记第i个区间内数据点的个数为f_i，令所有f_i除以n得到数据落入第i个区间的频率p_n，也就是

p_i = f_i/n

当n无穷大的时候，就可以认为p_i是数据在这个区间内的概率，并且有p_n=1，也就是具有归一性。

那么我们就直接做出区间右值与p_i的函数图象，就是这个数据的经验分布函数了。

写程序

假设我们现在的数据都在data里，是1xn大小的向量。

[~, n] = size(data);
%最大值与最小值
d_max = max(data);
d_min = min(data);
%划区间,step是小的区间长度
x = d_min:step:d_max;
len = length(x);
%提前给区间的频数f预留空间
func = zeros(1, len)
for i = 1:len
    %sum内的data<=(d_min+i*step)是逻辑判断语句
    %整个语句意思是找出data内小于等于d_min+i*step的总个数
    func(i) = sum(data<=(d_min+i*step))
end
plot(x, func/n, 'b-', 'LineWidth', 1)
title('经验分布函数')

第二步，做正态分布的图

要知道我们的经验函数虽然假设服从正态分布，所以我们要知道服从的是均值为多少，标准差为多少的正态分布，在这里我们就需要计算我们原始数据的均值mu与标准差sigma。

计算出mu和sigma后做出对应的正态分布图即可。

下面是程序

%均值
mu = sum(data)/n;
%标准差
sigma = sqrt(1/(n-1)*sum((data-mu).^2));
%正太分布数据点
func_norm = normcdf(x, mu, sigma);
plot(x, func_norm, 'k-', 'LineWidth', 1)

第三步，做置信带

实际上我也不是很清楚置信带具体的含义是什么，因为对于标准正态分布的95%置信带，意思就是样本点落在置信带内的概率为95%，但是如果是经验分布函数的话，画上一个置信带的意义就不是很明确了。

不过即便如此，这里还是画上比较好。

%f实际上就是经验分布函数，与上文中计算出来的func一模一样，x是横轴
%flo是置信下区间，fup是置信上区间
[f, x, flow, fup] = ecdf(data);
stairs(x, flow, 'r:', 'LineWidth', 1)
stairs(x, fup, 'r:', 'LineWidth', 1)

%写一下图例等
legend('经验分布函数', '正态分布函数', '置信下限', '置信上限')
xlabel('x')
ylabel('F(x)')

最后得出图像

如下图所示

可以看出来我们的经验分布函数和正态分布吻合的比较好。

感想

实际上我感觉置信带应该画正态分布的，而不是如果我们的经验分布函数所有点都在正态分布的置信带里面，则说明我们的数据服从正态分布，置信程度为95%。这样子还比较能够说得过去。