研读一本书,《数据之魅:基于开源工具的数据分析》(Data Analysis with Open Source Tools),写的很好。这里,复述一下书中用Gnuplot分析时间序列数据的部分。
Gnuplot安装很简单,直接到官网下载exe文件,安装运行即可(我是win7 32bit)。数据集来自这里,下载解压后,找到carbon.dioxide.txt文件,把内容复制到test.txt(方便而已),将test.txt文件放到Gnuplot的默认文件目录中(我的是在C:\Users\Administrator\Documents)。数据格式和书里面不太一样,只有一列,不过不影响分析。
首先,看下数据的总体趋势。
plot ‘test.txt‘ u 1 w l
可以解释为,画出,‘test.txt‘的图像,用(u-->using)第 1 列的数据,用直线表示(w l --> with line).
输出图形:
可以看出总体非线性的上升趋势和短期的周期性。
为了模拟方便,将数据减去一个垂直偏移,是图像经过原点。
plot ‘test.txt‘ u 0:($1-315) w l
这里,$1表是引用用数据的第一列的数据。用行号作为x坐标(位于伪列0),第一列(数据列)减去315,作为对应的y的坐标.
输出如下:
可以看到,曲线符合幂函数想x^k的特征,且曲线下凸,得出其k>1.
为了固定图的右上角(这里不太懂。。。),我们需要调整两轴,因为x^k通过(1,1),那么b(x/a)^k将通过(a,b)。这样就找到一个待定的拟合函数。点(35,350)在其附近,故得到函数35*(x/350)^k,这里仅知道k>1,经过测试,k=1.35时,拟合得最好。
plot ‘test.txt‘ u 0:($1-315) w l, 35*(x/350)**1.35
拟合图如下:
这样,我们继续,如果上面的函数拟合的足够好,那么消除这个总体趋势后,剩余部分将不会表现自身的趋势,将仅仅关于y=0对称。
plot ‘test.txt‘ u 0:($1-315-35*($0/350)**1.35) w l
从下面输出可以看出,基本是符合的。
在数据中很难看到长期趋势,一次我们用一天平滑曲线去逼近它。这里利用Gnuplot里面的加权样条法。
f(x) = 315 + 35*(x/350)**1.35 plot ‘test.txt‘ u 0:($1-f($0)) w l, "" u 0:($1-f($0)):(0.001) s acs w l
输出如下,这里的“0.001”越小,生成的曲线越平滑。
可以看到,整体剩余部分几乎是平坦的,横跨0的上下两端,也就是说,剩余部分几乎不存在系统偏差的迹象了。
完成对趋势的处理后,接下来要处理的特征就是周期性。
观察上图,可以看出其接近于正弦曲线,且振幅接近3.由于数据是按照月份统计的,所以我们猜测周期是一年。这意味着,每隔12个数据点,数据就是相同的。
plot ‘test.txt‘ u 0:($1-f($0)) w l, 3*sin(2*pi*x/12) w l
可以看到,数据拟合得特别好!同样的,让我们看看,本次拟合的剩余部分。
f(x) = 315 + 35*(x/350)**1.35 + 3*sin(2*pi*x/12) plot ‘test.txt‘ u 0:($1-f($0)) w l, "" u 0:($1-f($0)):(0.001) s acs w l
输出如下:
让我们放大一个具体的区间看一下。这里lp为线点类型,[60,120]自然是区间。
plot [60:120] ‘test.txt‘ u 0:($1-f($0)) w lp, "" u 0:($1-f($0)):(0.001) s acs w l
输出:
现在数出来两个主波谷间的数据点:12个!这是主要的季节性。但是,在两个主波谷之间还存在着一个次波谷,所以,图中含有更高次谐波(参考wiki解释,需FQ)。
图中看到,上升部分为7个月,下降部分为5个月,对于这种不对称性,我们试一下第一次高次谐波。
f(x) = 315 + 35*(x/350)**1.35 + 3*sin(2*pi*x/12) - 0.75*sin(2*pi*$0/6) plot ‘test.txt‘ u 0:($1-f($0)) w l, "" u 0:($1-f($0)):(0.001) s acs w l
输出:
现在,我们再来看下残差。
plot ‘test.txt‘ u 0:($1-f($0)) w l, "" u 0:($1-f($0)):(0.001) s acs w l,0,1,-1
输出:
观察到,残差偏向于正值,我们以0.1为单位调整垂直偏移。
f(x) = 315 + 35*(x/350)**1.35 + 3*sin(2*pi*x/12) - 0.75*sin(2*pi*$0/6) + 0.1 plot ‘test.txt‘ u 0:($1-f($0)) w l, "" u 0:($1-f($0)):(0.001) s acs w l,0,1,-1
输出:
可以看到,残差已经很小了。现在,我们结合分析模型再来看一下原始数据。
f(x) = 315 + 35*(x/350)**1.35 + 3*sin(2*pi*x/12) - 0.75*sin(2*pi*$0/6) + 0.1 plot ‘test.txt‘ u 0:1 w l , f(x)
输出:
Perfect!至此,我们已经得到一个解析式可以很好地描述数据。通过它,我们就可以相对准确地预测以后的数据。
plot [0:600] ‘test.txt‘ u 0:1 w l , f(x)
书中有段话写的很好:
那么问题的关键在哪里呢?关键在于,我们开始时对数据一无所知,即我们根本不知道数据会是什么样子。然后,经过逐层处理,我们剔除了数据的各个组成部分,直到只剩下随机噪音。最后,我们得到了一个明确的解析式,它很好地描述了数据。