高维空间随机统计软件HDS说明书.docx

高维空间随机统计软件HDS说明书北京精计软件科技有限公司第一章背景简介随着现代数据的收集和储存技术的提高，统计数据呈现出高维性。由于可重复研究的限制，参加研究的个体数量相对很小。这就是现代统计学中最具挑战的大P，小n问题。具体地说，数据的维数大大超过样本的个数。这尤其表现在生物基因学研究，网络信息，以及金融数据中。如何在样本量不是很大的前提下分析超高维数据，是一个非常具有挑战的、也是国际统计学的前沿课题。目前研究的重点从五个方面对高维数据统计建模与分析进行科学的，系统的研究。这五个方面是：（1）高维数据的变量选择、（2）超高维多元统计分析、（3）复杂数据的相关性、（4）大规模在线数据的监控和（5）高维生存数据分析。这五方面的研究均对传统的统计推断理论提出了全新的挑战，且均是目前国际统计学研究的最前沿问题。这五个课题相对独立又相互依托，有理论也有应用，将从不同的方向对高维数据的统计推断提出有效的解决方法，建立一个统一的适应于高维数据统计建模与分析的框架。本软件是一个高维空间随机统计软件，属于上述的第二个方向。目前已有的主要算法是马尔科夫链-蒙特卡洛（MCMC）方法，它可以搜索高维空间中的平衡分布，及其极值点。后面还将增加更多的快速、高效算法。对于MCMC算法，采样方法有多种，包括Metro-HaSting,GibbS等。版本LO为串行软件。初始参数点随机产生，一般经过一段时间的迭代后找到分布的近似极大值。MeMC是由两个MC构成的，分别指马尔科夫链和蒙特卡罗方法。马尔可夫链蒙特卡洛方法（MarkOVChainMonteCarlo）,简称MCMCo其产生于20世纪50年代早期，是在贝叶斯理论框架下，通过计算机进行模拟的蒙特卡洛方法（MOnteCarlo）o该方法将马尔可夫（MarkoV）过程引入到MonteCarIO模拟中，实现抽样分布随模拟的进行而改变的动态模拟，弥补了传统的蒙特卡罗积分只能静态模拟的缺陷。Metropolis等人在1953年首次提出了基于马尔可夫链的蒙特卡罗方法，即Metropolis算法，并在最早的计算机上编程实现。Metropolis算法是首个普适的采样方法，并启发了一系列MCMC方法，所以人们把它视为随机模拟技术腾飞的起点。MetroPOIiS算法这篇论文1被收录在统计学中的重大突破中，ComputinginScience&Engineering尝试列出了对20世纪科学与工程的发展和实践影响最大的十种算法。MetropolisAlgorithmforMonteCarlo被列为十大算法之首。用于蒙特卡洛的Metropolis算法定义了一个收敛的马尔可夫链，其极限就是所需的概率分布。Metropolis算法及其推广算法已被称为蒙特卡洛马尔可夫链技术(MCMC),因为这些算法模拟了一个马尔可夫链，从极限分布中获取抽样。1 .蒙特卡罗方法：是一种基于采样的随机近似方法，主要应用于随机采样、数学期望估计、定积分计算中。就是假设概率分布是已知的，然后通过采样获得概率分布的随机样本，得到了符合该概率分布的这些样本后，可以用于估计总体分布、计算均值来估计总体期望、通过期望计算积分等。所以蒙特卡罗方法核心就是随机采样。一般的蒙特卡罗方法有概率密度采样、接受拒绝采样。重要性采样等。2 .马尔科夫链：满足P(XtlXl-1,Xt-2,.,)=P(XtIXE),t=l,2.的随机序列乂=之”1,.”,.称为马尔科夫链，就是说XiE时的状态，只与Xt有关,而与之前的状态没有关系。马氏链在Xt“时的状态分布(t÷l)可以X时的状态分布(t)和转移概率矩阵P来决定，即(t+l)=P(t)o马氏链的一个重要性质就是平稳分布：若某时刻的一个状态分布使得JTP=JI(细致平稳方程()p(,*)=(x*)p(x*,x)是平稳分布的充分不必要条件)，则称兀为马尔可夫链的平稳分布。第二章算法原理1 .下面是最早的MH算法定义：1)初始化马氏链：Xo=Xo;2)对t=0,1,2,.进行下面的循环 第t个时刻马氏链状态为Xt=Xt,采样y=q(xxt)； 从均匀分布采样u=uniform(0.0,1.0); 如果u<(t,y)=p(y)q(xty),则接受转移Xt-y,即X+1=y; 否则不接受，即Xt+l=冗； 不断修正P(X)2 .下面是改进后的MH算法定义：1)初始化马氏链：X。=Xo;2)对t=0,1,2,进行下面的循环 第t个时刻马氏链状态为Xt=xt,采样y=q(t); 从均匀分布采样U=uniform(0.0,1.0); 如果U<a(xt,y)=minp(y)q(ly),1,则接受转移Xtfy,Pg)q(y%D即Xm=y; 否则不接受，即Xt+1=X1; 不断修正P(X)；3 .GibbS采样算法:1)初始化马氏链：Xo=%o；2)对t=0,l,2,进行下面的循环 %，+1)口(%11亚乩.,篇)； *+l)p(%2注+1，后，,篇)； 婿+l)p(%j悟+1,.,xjtt%j+1,.,Xn); X，l)P(Xn+1，厩+1，.,4D4.相似度函数形式：在调试MCMC算法时，采用的任意采样点信号同实际信号的相似度函数定义如下：122/1Slmmax(.jPOW(XLX?,2.0),1.0e-32)在GibbS算法中，由于只在一个方向进行采样，任意采样点信号同实际信号的相似度函数定义如下：sim=-器-(2)max(pow(Xj-X,2.0),1.0e-32)第三章软件结构1 .整体架构：下图是HDS软件的流程图。其中采样模块包含了样本点的采集和比对，以及根据算法的取舍。图LHDS的程序结构示意图2 .数据结构：下图是HDS采用的数据结构示意图。HDS软件图2.HDS的数据结构示意图(1) Sample数据类这是存储一个高维空间的样本点的类，包含了编号和空间坐标信息OClassSamplepublic:Sample(intdimension)coord=newintdimension;)Sampledntdimension,double*cds,double*bmin,double*dis);SampleOdeletecoord；Intid,ent;int*coord;；(2) SamPleSPaCe数据类这是存储所有高维空间的样本点的类，是高维随机统计的空间。其中有多个VeCtor和InaP类型的变量，用于进行样本点的采样、存储和排序等操作。它是进行高维随机统计的基础类。ClassSampleSpacepublic:SampleSpace(intns,intdimension,double*bmin,double*bmax,double*range,double*dis,int*mesh,double*target,double*start,intstype,intOutsampleO)dim=dimension;nsample=ns;coord=newD2D(nsample,dim);bmin=bmin;bmax=bmax;maxIoc=newdoubledim;mset=newstd:rmap<vector<int>,double>pmset=newMlD(dim);pmit=newstd:map<vector<int>,double>:iteratordim;e(nsample);；SampleSpace()deleteD2D(coord);deletemaxloc;free(mset);deleteMID(dim,pmset);deletepmit;；intnsample,tsample,dim;intoutsample;double*bmin,*bmax;double*maxloc,maxden;std:map<vector<int>,double>*mset;std:map<vector<int>,double>:iteratormit,mit;std:rmap<vector<int>,double>*pmset;std:map<vector<int>,double>:iterator*pmit;std:vector<int>vset;std:vector<int>:iteratorvit,wit;voidclearsample();voidcomputedensity();boolsearchsample(double*cds);voidinsertsample(double*cds);；(3) PSPaCe数据类这是并行环境下所有处理器上进行高维随机统计的类，SampleSPaCe是它的成员函数。其中进行样本点的随机搜集、随机采样、搜索和排序等算法。后面新的算法将增加更多的函数。classPSpacepublic:PSpace()rmax=1.0;rmin=0.0;sigma=0.2;gslrngenvsetup();T=gslrngdefault;r=gslrngalloc(T);sp=newSignal(r,10.0,1000.0);；PSpace()deletexmax;deleteparameters;deletebmin;deletebmax;deleterange;deletedis;deletemesh;deletetarget;deletestart;gslrngfree(r);deletetemp;deletefunc;free(gwp);；doublesigma,rmax,rmin;constgslrngtype*T;gslrng*r;intniter,nm;SampleSpace*data;GWtemplate*gwp;Signal*sp;voidreadin();voidinit();voidonesample(double*func);voidonesamplegibbs(intdir,double*func);doublesim(double*txyz);doublesim(intdir,double*txyz);doubleonedirgibbs(intdir,doubletsim,double*tden);doublerandstep(intdir);第四章软件结构和测试一.文件说明：1.主程序：1）头文件：mcmc.h：主程序头文件；util,h：工具函数头文件;2）源程序：mc2.C：主程序文件；sspace.C：样本空间类定义文件；pspace.C：并行空间类定义文件；input.C：参数读入函数文件；signal.C：信号处理函数定义文件；util.C：工具函数文件；random.C：随机函数测试文件；2 .调用了gsl库；3 ,调用了引力波模板库：FFT/Miscellaneous/ParametersMap/Wave/Wblock/PN/子目录下的文件都是引力波信号相关的文件。二.测试算例：1 .参数输入文件：（以8维为例如下）dimension8pari51.75e-l6.5e-l2.Oe-I0.0par251.75e-l6.5e-l-6.0e-20.0par351.75e-l6.5e-l1.Oe-I0.0par4 5 1. 75e-l 6. 5e-l1. 3e-l O. Opar5 5 1. 75e-l 6. 5e-l1. 5e-l O. Opar6 51. 75e-l 6. 5e-l -3. 5e-2 O. Opar7 5 1. 75e-l 6. 5e-l1. 7e-l O. Opar8 5 1. 75e-l6. 5e-l 2. 2e-l O. Osampling O/ O-Metropolis-Hasting;I-Gibbsnsample30000iteration300001outfileOOUtscreenOOUtsampleOnmodes16384/其中：dimension是维度数，下面跟着dimension行；每行包含6歹I。分别是维度名称（不重要）、这一维度上的网格数、下边界、上边界、测试极大值位置、初始位置；SamPIing可选O或1,O是指Metropolis-HaSting算法，1是指GibbS算法；nsample是指计算使用的样本数目iteration是指计算迭代的次数；OUtfile是指是否输出到文件；Outscreen是指是否输出到屏幕；OUtSanIPIe是指是否输出样本信息；nmodes是指信号波的采样点数；2. 三维测试：1）参数空间=0.175,0.653;2）参数空间网格：10（）3；3）测试参数点（参数空间中任意点）：0.2,-0.06,0.1；4）样本数二IOoO0;5）迭代次数=20000；6）计算结果：下图是3维空间两种采样方法的收敛历史测试。可以看到Gibbs采样方法比MH采样方法的效果更好，波动更小。ItefJtbn number图3.三维空间收敛测试3. 5维测试：1）参数空间二0.175,0.655;2）参数空间网格：205；3）测试参数点（参数空间中任意点）：0.2,-0.06,0.1,0.13,0.15；4）样本数=20000;5）迭代次数=IOooO0；6）计算结果：下图是5维空间两种采样方法的收敛历史测试。同样可以看到Gibbs采样方法比MH采样方法的效果更好，波动更小。由于维度增加，每一维度上的网格数减小，导致最终稿的误差增大。图4.5维空间收敛测4. 8维测试：1）参数空间=0.175,0.658;2）参数空间网格：58；3）测试参数点（参数空间中任意点）：0.2,-0.06,0.1,0.13,0.15,-0.035,0.17,0.42；4）样本数=30000;5）迭代次数=150000；6）计算结果：下图是8维空间两种采样方法的收敛历史测试。同样可以看到GibbS采样方法比MH采样方法的效果更好，波动更小。由于维度增加，每一维度上的网格数减小，导致最终的误差增大。而且随着维度的增加，样本数目和迭代次数明显增大，导致计算时间显著增大。0.1=F = = =*N=. jfc sfc 工尸 ? C Wf W 手 =25010 ' " 4l0 ' ' 75010 ' ' IO(M)IO ' ' ,2W1 , ' 1010IUfitIon ftps图5.8维空间收敛测试