计算数据的经验分布函数与MATLAB作图
写在前面
因为某些原因,需要处理某些数据,比如说某项测量数据与理论值的误差,我们就需要检验误差是否符合正态分布。最直观的方法就是直接做出经验分布函数的图来进行观察。所以这里简单的写一下如何做出经验分布图。
第一步,做经验分布图
对数据按升序排列
首先我们将所有一共n个数据按照升序排列。之所以升序,是因为这样方便我们后续计算。
将区间[a, b]分成小段
a是所有数据的最小值,b是所有数据的最大值,我们将所有数据分成小区间,比如说将[a,b]分为
- [a,a+(b-a)/100]
- [a, a+2*(b-a)/100]
- ...
- [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%。这样子还比较能够说得过去。