matlab__数据分析和统计资料

1、Good is good, but better carries it.精益求精，善益求善。matlab_数据分析和统计-数据分析和统计面向列的数据集这年头似乎十分风行”面向”这个词，这儿故也套用，其英文为Column-OrientedDataSets,可理解为MatLab按列的存储方式来分析数据,下面是一个例子:TimeLocation1Location2Location301h001111902h007131103h0014172004h001113905h0043516906h0038467607h006113218608h007513518009h00388811510h0028365

2、511h0012121412h0018273013h0018192914h0017151815h0019364816h0032471017h0042659218h00576615119h0044559020h0011414525721h0035586822h0011121523h001391524h001097以上数据被保存在一个称为count.dat的文件中11119713111417201113943516938467661132186751351803888115283655121214182730181929171518193648324710426592576615144559011

3、4145257355868111215139151097下面，我们调入此文件，并看看文件的一些参数loadcount.datn,p=size(count)n=24p=3创建一个时间轴后，我们可以把图画出来：t=1:n;set(0,defaultaxeslinestyleorder,-|-|-.)set(0,defaultaxescolororder,000)plot(t,count),legend(Location1,Location2,Location3,0)xlabel(Time),ylabel(VehicleCount),gridon足以证明，以上是对个对象的次观测基本数据分析函数（一定

4、注意是面向列的）继续用上面的数据，其每列最大值均值及偏差分别为：mx=max(count)mu=mean(count)sigma=std(count)mx=114145257mu=32.000046.541765.5833sigma=25.370341.405768.0281重载函数，还可以定位出最大最小值的位置mx,indx=min(count)mx=797indx=22324试试看，你能看懂下面的命令是干什么的吗？n,p=size(count)e=ones(n,1)x=counte*muHYPERLINKjavascript:void(null)点这看看答案！下面这句命令则找出了整个矩阵的

5、最小值:min(count(:)ans=7协方差及相关系数下面,我们来看看第一列的方差:cov(count(:,1)ans=643.6522cov()函数作用于矩阵,则会计算其协方差矩阵.corrcoef()用于计算相关系数,如:corrcoef(count)ans=1.00000.93310.95990.93311.00000.95530.95990.95531.0000数据的预处理未知数据NaN(NotaNumber-不是一个数)被定义为未经定义的算式的结果,如0/0.在处理数据中,NaN常用来表示未知数据或未能获得的数据.所有与NaN有关的运算其结果都是NaN.a=magic(3);a(

6、2,2)=NaNa=8163NaN7492sum(a)ans=15NaN15在做统计时，常需要将NaN转化为可计算的数字或去掉,以下是几种方法:注:判断一个值是否为NaN,只能用isnan(),而不可用x=NaN;i=find(isnan(x);x=x(i)先找出值不是NaN的项的下标,将这些元素保留x=x(find(isnan(x)同上,去掉NaNx=x(isnan(x);更快的做法x(isnan(x)=;消掉NaNX(any(isnan(X),:)=;把含有NaN的行都去掉用此法可以从数据中去掉不相关的数据,看看下面的命令是干什么用的:mu=mean(count);sigma=std(co

7、unt);n,p=size(count)outliers=abs(countmu(ones(n,1),:)3*sigma(ones(n,1),:);nout=sum(outliers)nout=100count(any(outliers),:)=;HYPERLINKjavascript:void(null)点这看看答案回归与曲线拟合我们经常需要把观测到的数据表达为函数，假如有如下的对时间的观测：t=1.62.3;y=0.50.851.40;plot(t,y,o),gridon多项式回归由图可以看出应该可以用多项式来表达：y=a0+a1*t+a2*t2系数a0,a1,a2可以由最小平方拟合来确定

8、，这一步可由反除号来完成解下面的三元方程组可得：X=ones(size(t)tt.2X=1.0000001.00000.30000.09001.00000.80000.64001.00001.10001.21001.00001.60002.56001.00002.30005.2900a=Xya=0.53180.91910.2387a即为待求的系数，画图比较可得T=(0:0.1:2.5);Y=ones(size(T)TT.2*a;plot(T,Y,t,y,o,),gridon结果令人失望，但我们可以增加阶数来提高精确度，但更明智的选择是用别的方法线性参数回归形如：y=a0+a1*exp(-t)+

9、a2*t*exp(-t)计算方法同上：X=ones(size(t)exp(t)t.*exp(t);a=Xya=1.39740.89880.4097T=(0:0.1:2.5);Y=ones(size(T)exp(T)T.exp(T)*a;plot(T,Y,t,y,o),gridon看起来是不是好多了！例子研究：曲线拟合下面我们以美国人口普查的数据来研究一下有关曲线拟合的问题（MatLab是别人的,教学文档是别人的,例子也是别人的，我只是一个翻译而已）loadcensus这样我们得到了两个变量，cdate是1790至1990年的时间列向量(10年一次),pop是相应人口数列向量.上一小节所讲的多项

10、式拟合可以用函数polyfit()来完成,数字指明了阶数p=polyfit(cdate,pop,4)Warning:Matrixisclosetosingularorbadlyscaled.Resultsmaybeinaccurate.RCOND=5.429790e20p=1.0e+05*0.00000.00000.00000.01266.0020产生警告的原因是计算中的cdata值太大，在计算中的Vandermonde行列式使变换产生了问题，解决的方法之一是使数据标准化预处理：标准化数据数据的标准化是对数据进行缩放，以使以后的计算能更加精确，一种方法是使之成为0均值：sdate=(cdate

11、mean(cdate)./std(cdate)现在再进行曲线拟合就没事了！p=polyfit(sdate,pop,4)p=0.70470.921023.470673.859862.2285pop4=polyval(p,sdate);plot(cdate,pop4,cdate,pop,+),gridon在上面的数据标准化中，也可以有别的算法，如令年的人口数为余量分析p1=polyfit(sdate,pop,1);pop1=polyval(p1,sdate);plot(cdate,pop1,cdate,pop,+)res1=poppop1;figure,plot(cdate,res1,+)p=po

12、lyfit(sdate,pop,2);pop2=polyval(p,sdate);plot(cdate,pop2,cdate,pop,+)res2=poppop2;figure,plot(cdate,res2,+)p=polyfit(sdate,pop,4);pop4=polyval(p,sdate);plot(cdate,pop4,cdate,pop,+)res4=poppop4;figure,plot(cdate,res4,+)可以看出，多项式拟合即使提高了阶次也无法达到令人满意的结果指数拟合从人口增长图可以发现人数的增长基本是呈指数增加的，因此我们可以用年份的对数来进行拟合，这儿，年数是

13、标准化后的！logp1=polyfit(sdate,log10(pop),1);logpred1=10.polyval(logp1,sdate);semilogy(cdate,logpred1,cdate,pop,+);gridonlogres2=log10(pop)polyval(logp2,sdate);plot(cdate,logres2,+)上面的图不令人满意，下面，我们用二阶的对数分析：logp2=polyfit(sdate,log10(pop),2);logpred2=10.polyval(logp2,sdate);semilogy(cdate,logpred2,cdate,pop

14、,+);gridonr=pop10.(polyval(logp2,sdate);plot(cdate,r,+)这种余量分析比多项式拟合的余量分析图案要随机的多（没有很强的规律性），可以预见，随着人数的增加，余粮所反映的不确定度也在增加，但总的说来，这种拟合方式要强好多！误差边界误差边界常用来反映你所用的拟合方式是否适用于数据，为得到误差边界，只需在polyfit()中传递第二个参数，并将其送入polyval().下面是一个二阶多项式拟合模型,年份已被标准化,下面的代码用了2,对应于95%的可置信度:p2,S2=polyfit(sdate,pop,2);pop2,del2=polyval(p2,

15、sdate,S2);plot(cdate,pop,+,cdate,pop2,g,cdate,pop2+2*del2,r:,.cdate,pop22*del2,r:),gridon差分方程和滤波MatLab中的差分和滤波基本都是对向量而言的，向量则是存储取样信号或序列的函数y=filter(b,a,x)将用a,b描述的滤波器处理向量x,然后将其存储在向量y中,filter()函数可看为是一差分方程a1y(n)=b1*x(1)+b2*x(2)+.-a2*y(2)-.如有以下差分方程:y(n)=1/4*x(n)+1/4*x(n-1)+1/4*x(n-2)+1/4*x(n-3),则a=1;b=1/41

16、/41/41/4;我们载入数据,取其第一列,并计算有:loadcount.datx=count(:,1);y=filter(b,a,x);t=1:length(x);plot(t,x,.,t,y,),gridonlegend(OriginalData,SmoothedData,2)实现所表示的就是滤波后的数据,它代表了4小时的平均车流量MatLab的信号处理工具箱中提供了很多用来滤波的函数,可用来处理实际问题!快速傅立叶变换(FFT)傅立叶变换能把信号按正弦展开成不同的频率值,对于取样信号,用的是离散傅立叶变换.FFT是离散傅立叶变换的一种高速算法,在信号和图像处理中有极大的用处!fft离散傅

17、立叶变换fft2二维离散傅立叶变换fftnn维离散傅立叶变换ifft离散傅立叶反变换ifft2二维离散傅立叶反变换ifftnn维离散傅立叶反变换abs幅度angle相角unwrap相位按弧度展开,大于的变换为2的补角fftshift把零队列移至功率谱中央cplxpair把数据排成复数对nextpow2下两个更高的功率向量x的FFT可以这样求:x=43791000y=fft(x)y=6.000011.48532.7574i2.000012.0000i5.4853+11.2426i18.00005.485311.2426i2.0000+12.0000i11.4853+2.7574ix虽然是实数,但

18、y是复数,其中,第一个是因为它是常数相加的结果,第五个则对应于奈奎斯特频率,后三个数是由于负频率的影响,它们是前面三个数的共轭值!下面,让我们来验证一下太阳黑子活动周期是11年!Wolfer数记录了300年太阳黑子的数量及大小:loadsunspot.datyear=sunspot(:,1);wolfer=sunspot(:,2);plot(year,wolfer)title(SunspotData)现在来看看其FFT:Y=fft(wolfer);Y的幅度是功率谱,画出功率谱和频率的对应关系就得出了周期图,去掉第一点,因为他只是所有数据的和,画图有:N=length(Y);Y(1)=;power=abs(Y(1:N/2).2;nyquist=1/2;freq=(1:N/2)/(N/2)*nyquist;plot(freq,power),gridonxlabel(cycles/year)title(Periodogram)上面的图看起来不大方便,下面我们画出频谱-周期图period=1./freq;plot(period,power),axis(04002e7),gridonylabel(Power)xlab