GnuPlot：简单数据统计处理及取整函数int/floor的问题
2012-11-17 09:27:14 初稿
2013-08-02 14:08:45 修订
2014-04-08 20:02:45 再订

最近需要利用GnuPlot对数据进行简单的统计处理，就查阅了一些相关资料，在此简单总结一下我所用到的方法。

本质上说，GnuPlot只是个画图软件，画图功能是其强项，而计算、统计本不属于它的功能，但是利用GnuPlot的函数及脚本功能，我们还是可以变通地利用它进行一些计算及统计处理，这些都属于tricks，若感兴趣，请参看此处的总结，取其所需。

上面所引网页对数据的统计处理是利用GnuPlot的内置变量，略显罗嗦，不太合我的胃口。我喜欢利用自定义函数来处理，方式如下：

●    数据的最值：定义两个变量min和max，再定义一个函数minmax
min=1E38; max=-1E38;
minmax(x)=(x<min)?min=x:(x>max)?max=x:0
plot "file" u 1:(minmax($1)); print min, max

●    数据的均值，标准差：定义变量avg和函数mean
mean(x)=avg
fit mean(x) "file" u 1:1 via avg

根据fit的输出直接就可以得到平均值avg和标准差.

新版本的GnuPlot, 可以直接使用stat命令得到数据的统计信息, 不需要上面这么麻烦.

当然，如果你只是做一次两次，像我这么做并没有什么优势，但是如果你要经常使用，那最好还是把这些函数添加到GnuPlot的脚本中，这样就可以直接使用了，就像那些内置函数一样。

如何利用GnuPlot做统计直方图（Histogram），这是经常有人问的问题。对此，官方有答案，stackoverflow有答案，马欢也简单说过。实现方法很简单，定义一个binning函数，bin(x,w)=w*int(x/w)或bin(x,w)=w*floor(x/w)（若希望画图时数据处于区间中心，可使用bin(x,w)=w*(int(x/w)+0.5)或bin(x,w)=w*(floor(x/w)+0.5)），w为bin间隔大小。此函数可以把数据映射到相应的区间，再利用GnuPlot的smooth功能进行频数统计，画出来即可。具体来说就是

概率密度：Plot “DataFile” u(bin($1,w)):(1./Ndat/w) s f
累积分布：Plot “DataFile” u 1:(1./Ndat) s cum

Ndat为数据个数，w为bin的宽度，1./Ndat/w用于归一化。

有时，我们已经有了某个分布的概率密度数据，但为了作图平滑，我们希望改变bin的宽度，这种情况下Ndat是不变的，改变的只是w，可使用

改变bin宽度：Plot “DataFile” u (bin(\$1,w)):(\$2*w0/w)s f

w0为原先的bin宽度，w为新的宽度。

本来是很简单的事情，但我使用时却发现了问题，统计出来的频数有时候不正常。究其原因，在于int和floor这两个内置的取整函数在某些情况下给出的结果和精确结果差1。如int(0.0003/0.0001)=floor(0.0003/0.0001)=2。在数据量很大的情况下，这点差距影响并不大，但在数据量比较小的情况下，很可能对结果产生较大影响。我进而测试了awk和perl的int函数，给出了同样的结果。这使我认识到这个问题有其内在原因。查一下，果然，请看这里。原来是因为计算机存储数据时对某些数据并不能精确存储，如计算0.0003/0.0001，得到的结果并不是精确的3，很可能是2.99999999999999955591，这样用int或floor取整后就变成2了。

既然知道了原因，那解决也很简单，有人提到定义一个递归bin函数来解决，我觉得过于麻烦，计算也慢，我们直接改造int函数就可以了。重新定义一个rint函数，rint(x)=(x-int(x)>0.9999)?int(x)+1:int(x)，利用这个函数进行取整就没问题了。

上面定义的rint函数使用时对负数存在问题, 可以利用nint函数解决. GnuPlot没有内置这个函数, 但可利用floor和ceil函数的组合来定义

nint(x)=(x>0.)?floor(x+0.5):ceil(x-0.5)

bin(x,s)=s*nint(x/s)

使用这个bin函数统计就不会有问题了.