第十章数据处理对于实际的数学问题,我们经常需要处理各种各样的数据,在这一章中,我们将介绍Maple 处理数据的各种方法.§10.1统计程序包Maple的统计程序包由7个子程序包和一个命令组成.它们分别是:•anova方差分析程序包.•describ描述统计程序包.•fit线形回归程序包.•random随机数生成程序包.•statevalf统计分布计算程序包.•statplots统计图形程序包.•transform数据处理程序包.•importdate导入数据文件.下面我们通过例子来说明stats的用法.首先输入一些数据>with(stats);[anova,describe,fit,importdata,random,statevalf,statplots,transform] >with(describe);[coefficientofvariation,count,countmissing,covariance,decile,geometricmean, harmonicmean,kurtosis,linearcorrelation,mean,meandeviation,median,mode,moment, percentile,quadraticmean,quantile,quartile,range,skewness,standarddeviation,sumdata,variance]>data:=[2,3.2,4,Weight(3,5),missing,Weight(1..4,3)];data:=[2,3.2,4,Weight(3,5),missing,Weight(1..4,3)]>mean(data);2.881818182在我们输入数据的过程中missing代表缺失的数据,Weight代表分量.其中Weight(3,5)表示分量3重复了5次,而Weight(1..4,3)表示1到4之间的3个数据,因此是(1+4)/2×3.在实际应用中,输入数据的方式有许多种,在Maple中,我们还可以用importdata从文件中导入数据.例如给定数据文件statdata.dat,它包含的数据为12 3.51*4273987.5229230第十章数据处理在使用importdata命令时如果不加参数,它将把数据按照一个序列读出.>importdata(‘c:\\statdata.dat‘);1.,2.,3.5,1.,missing,4.,2.,7.,3.,9.,8.,7.5如果加上一个参数,它将把数据按照指定的列数读出.>importdata(‘c:\\statdata.dat‘,3);[1.,1.,2.,9.],[2.,missing,7.,8.],[3.5,4.,3.,7.5]对于输入的数据,transform子程序包提供了许多命令进行处理>with(transform);[apply,classmark,cumulativefrequency,deletemissing,divideby,frequency,moving, multiapply,scaleweight,split,standardscore,statsort,statvalue,subtractfrom,tally,tallyinto]>tst:=[Weight(1..3,4),3..4,Weight(4..7,4),8];tst:=[Weight(1..3,4),3..4,Weight(4..7,4),8]statvalue命令可以把输入的数据中的Weight命令移除,并按照Weight的含义进行转换.使用classmark命令可以把所有的表示范围的数据换为它们的平均值.frequency命令把Weight类型的数据替换为它的频数.cumulativefrequency命令的作用和frequency相似,但它把频数进行了累加.>statvalue(tst);[1..3,3..4,4..7,8]>classmark(tst),statvalue(classmark(tst));[Weight(2,4),72,Weight(112,4),8],[2,72,112,8]>frequency(tst);[4,1,4,1]>cumulativefrequency(tst);[4,5,9,10]这里的apply作用和map的作用是相同的,但是Weight因子不变.multiapply的作用和zip 相同,也是生成一个新的列表.>apply[sin](tst);[Weight(sin(1)..sin(3),4),sin(3)..sin(4),Weight(sin(4)..sin(7),4),sin(8)]>multiapply[(x,y)->x*y]([[1,2,3],[4,5,6]]);[4,10,18]§10.1统计程序包231scaleweight命令把方括号内的数与Weight因子相乘,在下面的例子中,count是describe 子程序包中的命令,用来对数据列表记数.>scaleweight[1/count(tst)](tst);[Weight(1..3,25),Weight(3..4,110),Weight(4..7,25),Weight(8,110)]statsort可以把数据列表进行排序.>statsort([6,missing,3,Weight(4..5,2)]);[3,Weight(4..5,2),6,missing]split命令把数据列表按照方括号内的数字n分割为等长的n份,但是数据的次序不变.如果数据的个数不能n等分,split将在列表中插入一些数据(使用Weight因子).>split[3]([3,4,5,6,7,8]);[[3,4],[5,6],[7,8]]>split[4]([3,4,5,6,7,8]);[[3,Weight(4,12)],[Weight(4,12),5],[6,Weight(7,12)],[Weight(7,12),8]]tallyinto命令的作用是按照给定的模式对数据进行分类,并统计数据出现的频数.返回的结果是Weight类型的数据.在下面的例子中,第一项为Weight(1..3,4)是因为数据中1≤x<3的x共有4项.>tallyinto([1,2,2,2,3,4,5,6,6],[1..3,3..4,4..7]);[Weight(1..3,4),3..4,Weight(4..7,4)]>tst:=[seq(rand(100)(),n=1..100)];tst:=[81,70,97,63,76,38,85,68,21,9,55,63,57,60,74,85,16,61,7,49,86,98,66, 9,73,81,74,66,73,42,91,93,0,11,38,13,20,44,65,91,95,74,9,60,82,92,13,77, 49,35,61,48,3,23,95,73,89,37,57,99,94,28,15,55,7,51,62,97,88,42,97,98,27,27,74,25,7,82,29,52,4,85,45,98,38,76,75,74,23,0,19,1,49,47,13,65,44,11,36,59]>bins:=[seq(10*n..10*n+10,n=0..9)];bins:=[0..10,10..20,20..30,30..40,40..50,50..60,60..70,70..80,80..90,90..100]>tallyinto(tst,bins);[Weight(30..40,6),Weight(40..50,10),Weight(50..60,7),Weight(60..70,12),Weight(70..80,13),Weight(80..90,10),Weight(90..100,14),Weight(0..10,11),Weight(10..20,8),Weight(20..30,9)]232第十章数据处理>frequency(%);[6,10,7,12,13,10,14,11,8,9]在统计程序包中,describe子程序包提供了一些常规的统计计算命令.例如:count可以计算输入数据的个数,countmissing则能给出输入数据中缺失数据的个数.range可以求出数据的范围,mean给出平均值.>tst:=[1,2,3,Weight(4..6,3),Weight(missing,2)]:>count(tst);6>countmissing(tst);2>range(tst);1..6>mean(tst);72计算其它平均值的命令还有quadraticmean,harmonicmean和geometricmean.计算标准差和方差的命令是standarddeviation和variance,它们的变种还有coefficientofvariation 和meandeviation.>quadraticmean(tst);16√534>standarddeviation(tst);16√93>variance(tst);3112对于多组数据的比较,可以使用的命令是linearcorrelation和covariance. >tst:=[1,2,3,4],[2,4,5,6]:>linearcorrelation(tst);1335√7>evalf(%);.9827076297§10.1统计程序包233 >covariance(tst);evalf(%);1381.625000000对于其它描述统计的命令读者可以通过Maple的帮助学习,在这里就不做详细的说明了.下面我们介绍statplots子程序包的用法.statplots的命令主要用来画各种统计图形,我们首先看一下直方图.>with(stats):with(statplots):with(transform):>with(random):with(describe):with(plots):>f:=normald[5,1]:>data:=[f(100)]:为了画直方图,我们首先生成一组数据,命令normald[5,1]定义了一个满足正态分布的随机分布函数f,其中µ=5,σ=1.下一个命令使用函数f生成了100个随机数据.对于其它的随机分布,可以用类似的方式得到,有关细节可参照Maple的帮助.>tallyinto(data,[seq(n*0.5..(n+1)*0.5,n=1..20)]);[Weight(9.0..9.5,0),Weight(8.5..9.0,0),Weight(6.0..6.5,4),2.5..3.0,Weight(4.0..4.5,15),Weight(9.5..10.0,0),Weight(6.5..7.0,4),Weight(3.0..3.5,4),Weight(.5..1.0,0),Weight(10.0..10.5,0),Weight(7.0..7.5,4),Weight(4.5..5.0,21),Weight(3.5..4.0,14),Weight(1.0..1.5,0),7.5..8.0,Weight(5.0..5.5,12),Weight(1.5..2.0,0),Weight(8.0..8.5,0),Weight(5.5..6.0,20),Weight(2.0..2.5,0)] >histdata:=scaleweight[1.0/count(")](");histdata:=[Weight(9.0..9.5,0),Weight(8.5..9.0,0),Weight(6.0..6.5,.04000000000), Weight(2.5..3.0,.010********),Weight(4.0..4.5,.1500000000),Weight(9.5..10.0,0),Weight(6.5..7.0,.04000000000),Weight(3.0..3.5,.04000000000),Weight(.5..1.0,0),Weight(10.0..10.5,0),Weight(7.0..7.5,.04000000000),Weight(4.5..5.0,.2100000000),Weight(3.5..4.0,.1400000000),Weight(1.0..1.5,0),Weight(7.5..8.0,.010********),Weight(5.0..5.5,.1200000000),Weight(1.5..2.0,0),Weight(8.0..8.5,0),Weight(5.5..6.0,.2000000000),Weight(2.0..2.5,0)]使用tallyinfo和scaleweight命令对数据进行分类处理以后,就可以用histogram命令画出直方图.但是为了同时画出正态分布曲线,我们先用p1记录这个直方图.画正态分布曲线的方法是使用statevalf子程序包,这个子程序包的功能是计算各种统计分布函数的值.它的一般用法是statevalf[function,distribution](args).对于连续的分布,function可取的值及含义如下:cdf累积密度函数icdf逆累积密度函数pdf概率密度函数对于离散的分布,function可取的值及含义如下:234第十章数据处理dcdf 离散累积概率函数idcdf 逆离散累积概率函数pf 概率函数statevalf 的第二个参数是distribution ,它可以是任何一种统计分布函数,例如我们前面定义的分布函数f .下面的命令fpdf:=statevalf[pdf,f]定义了一个函数,用它可以画出µ=5,σ=1的正态分布曲线.>p1:=histogram(histdata):>fpdf:=statevalf[pdf,f]:>p2:=plot(fpdf(x),x=0.5..10):>display([p1,p2]);使用statplots 子程序包中的scatter2d 和quantile2命令可以画二维数据的散点图.点的坐标来自于两个列表.下面我们先建立两个列表,scatter2d 以这两个列表作为参数画出散点图.>datax:=[normald[5,1](100)]:>datay:=[normald[7,3](100)]:>scatter2d(datax,datay);0246810121416234567quantile2命令的用法和scatter2d 的用法相同,唯一的区别是,quantile2先把两个列表的数据进行从小到大的排序,然后在画出散点图.>quantile2(datax,datay);§10.1统计程序包2350246810121416234567对于一维数组,我们可以用scatter1d ,boxplot ,notchedbox 等命令来画图.>scatter1d(datax);00.511.52234567scatter1d 命令将所有的点画在y =1的轴上.当使用jittered 选项时,scatter1d 会在0到1之间随机选择y 的值来画散点图.>scatter1d[jittered](datax);00.20.40.60.8234567boxplot 命令可以用箱式图描述数据的分布情况,notchedbox 的作用和boxplot 相似,其方括号内的参数可以指定箱子的纵向位置和宽度.>boxplot(datax);236第十章数据处理234567-0.4-0.200.20.4>notchedbox[5,2](datax);2345674 4.55 5.5在statplots 子程序包中还有三个命令分别是:xshift ,xscale 和xyexchange 它们都可以作用到二维图形上.xshift 的作用是在x 轴方向平移图形,xscale 的作用是对x 值进行缩放,xyexchange 命令可以交换x 和y 的坐标,作用到二维图形上就是旋转90◦.下面我们画一个二维图形,然后对它进行平移、缩放、旋转.>P1:=plot(x^2,x=0..3):>P2:=xshift[3](P1):>P3:=xscale[2](P1):>P4:=xyexchange(P1):>plots[display]([P1,P2,P3,P4]);00.511.522.5312345x§10.2插值与回归237§10.2插值与回归插值的含义是给出一个光滑的曲线让它通过给定的点.在Maple 中有两个命令处理插值问题,它们分别是interp 和spline .interp 可以求出一个多项式,使它恰好通过指定的点.例如,我们首先随机给出10个点>datax:=[seq(i,i=1..10)];datax :=[1,2,3,4,5,6,7,8,9,10]>datay:=[seq(rand(10)(),i=1..10)];datay :=[1,0,7,3,6,8,5,8,1,9]>dataxy:=zip((x,y)->[x,y],datax,datay);dataxy :=[[1,1],[2,0],[3,7],[4,3],[5,6],[6,8],[7,5],[8,8],[9,1],[10,9]]使用interp 命令可以求出一个次数为9的多项式f>f:=interp(datax,datay,x);f :=1751840x 9−51740320x 8+1169960480x 7−37192880x 6+2732317280x 5+1767415760x 4−6525773240x 3+18164833360x 2−166********x +293然后我们画出这条曲线以及10个点>p1:=plot(f,x=0.9..10.1):>p2:=plot(dataxy,style=point,symbol=circle,color=blue):>plots[display]([p1,p2]);-10-55101520252468x 虽然interp 求出的函数非常光滑,但是它的次数太高,要用这个函数进行外推的话偏差太大.另外一种插值方法是样条插值,命令为spline ,它采用分段插值的方法来得到一个分段光滑的曲线,通常每一段都是2次或3次的曲线.下面我们就用上面给出的数据进行样条插值.>readlib(spline):由于spline 不是Maple 中的常用命令,所以我们需要用readlib 命令来调用它.spline 的前两个参数分别为x 和y 坐标的数据列表,第三个参数为函数的变量,第四个参数为分段曲238第十章数据处理线的次数,通常我们使用2次或3次的曲线,但是spline 也可以计算更高次数的曲线.下面我们给出了2次和3次样条的函数和图象.>g:=unapply(spline(datax,datay,x,2),x);g :=x →piecewise(x <2,2x −x 2,x <3,40−38x +9x 2,x <4,−221+136x −20x 2,x <5,531−240x +27x 2,x <6,−844+310x −28x 2,x <7,992−302x +23x 2,x <8,−968+258x −17x 2,x <9,568−126x +7x 2,649−144x +8x 2)>h:=unapply(spline(datax,datay,x,3),x);h :=x →piecewise(x <2,2+20780940545x −12417713515x 2+12417740545x 3,x <3,114890213515−9681238109x +71722713515x 2−593058109x 3,x <4,−8287872703+6492882385x −104705213515x 2+29156840545x 3,x <5,10123012703−967434440545x +67896813515x 2−139********x 3,x <6,3143262703−339256940545x +26018313515x 2−212153x 3,x <7,−1705167413515+2454238740545x −129175913515x 2+20247740545x 3,x <8,778912255−5045848340545x +227971113515x 2−30773340545x 3,x <9,−253485124505+8153537340545x −322003313515x 2+759478109x 3,4282833901−5875265840545x +3951642703x 2−19758240545x 3)>p3:=plot(g(x),x=1..10):>p4:=plot(h(x),x=1..10,color=green):>plots[display]([p2,p3,p4]);-2024681012142468x 回归的含义是构造一条近似的曲线,使它尽可能的接近给定的数据点.在统计程序包的fit 子程序包中的leastsqr 命令可以用来构造回归函数.下面我们通过例子来说明leastsqr 的用法.首先我们生成一组随机数据,然后画出它们的散点图p1.§10.2插值与回归239>with(stats):with(fit):with(statplots):with(plots):>datax:=[seq(rand(100)(),i=1..50)]:>datay:=map(x->evalf(sin(x/100*Pi)+rand(30)()/100,4),datax):>p1:=scatter2d(datax,datay):我们用一个二次曲线y =a 0+a 1∗x +a 2∗x 2来拟合这些数据.使用的命令是leastsquare .它的参数有两组.在方括号中的第一组参数用来描述拟合曲线,包括:变量列表,拟合函数形式,以及拟合函数的参数集合.在圆括号中的第二组参数是数据的列表.具体用法如下:>leastsquare[[x,y],y=a0+a1*x+a2*x^2,{a0,a1,a2}]([datax,datay]);y =.1866122788+.0391*******x −.0003984139193x 2>f:=rhs(%);f :=.1866122788+.0391*******x −.0003984139193x 2拟合函数f 的曲线用p2表示.用display 命令同时显示p1和p2.>p2:=plot(f,x=0..100):>display([p1,p2]);0.20.40.60.811.2020*********x leastsquare 还可以处理多维的简单拟合函数.例如y =a 1∗x 1+a 2∗x 2+···.>datax1:=[1,2,3]:datax2:=[2,3,4]:datay1:=[3,7,8]:>leastsquare[[x1,x2,y],y=a1*x1+a2*x2,{a1,a2}]([datax1,datax2,datay1]);y =32x1+x2需要说明的是:leastsquare 命令使用的是最小二乘法.它对于拟合函数有一些限制.首先它要求拟合的参数必须是线性的.对于非线性的函数,例如f =a 0+a 1∗sin(a 2∗x )leastsquare 命令无法使用.此外拟合函数本身也必须是比较简单的函数,通常是多项式函数.对于复杂的函数,例如:f =a 0+a 1∗sin(x/10)+a 2∗sin(x/20)+a 3∗sin(x/40),leastsquare 虽然也可以计算,但是经过很长时间以后也未必能返回正确的结果.这样的问题是从Maple V Release 4以后的版本才出现的.在Maple V Release 3版本中,还可以求出正确的结果.不过使用共享库中的fit 命令,还是可以解决上述问题的.下面我们就介绍一下fit 命令的用法.首先我们调入共享库.>with(share):>readshare(fit,numerics):240第十章数据处理对于第一个例子的数据,fit 命令的用法如下:>f1:=fit(datax,datay,[1,x,x^2],x);f1:=.1866122947+.0391*******x −.0003984139092x 2可以看出计算的结果和leastsquare 命令是一样的.对于复杂的非线性函数,fit 同样有效,例如我们的拟合函数为f =a 0+a 1∗sin(x/10)+a 2∗sin(x/20)+a 3∗sin(x/40).计算的结果如下:>f:=fit(datax,datay,a0+a1*sin(x/10)+a2*sin(x/20)+a3*sin(x/40),x);f :=.0321572783+.07881337051sin(110x )+.3489868374sin(120x )+1.007798304sin(140x )>p3:=plot(f,x=0..100):>display([p1,p3]);0.20.40.60.811.2020*********x fit 也可以处理多维的情况.例如:>datax1:=[seq(rand(100)()/100,i=1..50)]:>datax2:=[seq(rand(100)()/100,i=1..50)]:>datax12:=zip((x1,x2)->[x1,x2],datax1,datax2):>datay:=zip((x1,x2)->x1+2*x2+cos(x1*x2)+rand(100)()/1000,datax1,datax2):>f:=fit(datax12,datay,a1*x1+a2*x2+a3*cos(x1*x2)+a4*x1*x2,[x1,x2]);f :=1.058776261x1+2.055120765x2+1.000361369cos(x1x2)−.05052181835x1x2§10.3极大与极小Maple 可以使用maximize ,minimize 和extrema 命令计算一维和多维函数的极大值和极小值.这三个命令都需要用readlib 来调入.我们首先介绍maximize 和minimize 的使用方法.>readlib(minimize):readlib(maximize):>minimize((x+2)/(x^2+1));1−12√5§10.3极大与极小241>minimize(sin(x));−1对于多维函数minimize 命令同样可以使用.>minimize(x^2+y^2+5);5对于多维函数,我们还可以指定对某个变量求极值.>minimize(x^2+3*y^3-1,{x });3y 3−1此外还可以在一个指定的范围内求极值.>minimize(x^4+2*y^2*x,{x,y },{x=-10..10,y=2..5});752RootOf(2Z 3+25,−2.320794417)>evalf(%);−87.02979064在Maple V Release 4版本中,minimize 有一个选项’infinite’.对于那些无法求出极值的函数,使用这个选项会返回函数自身.例如函数e x 是没有极值的,因此minimize 不返回值,但是使用了’infinite’选项后,它会返回e x 自身.不过在Maple 6以后的版本中,minimize(exp(x)得到的结果是0,它实际上是lim x →−∞e x 的极限值.>minimize(exp(x));>minimize(exp(x),’infinite’);e x对于maximize 命令,它实际上是把表达式加了减号,然后取极小值,再对极小值取减号得到的结果.因此我们不再举例来说明它的用法.在Maple 6以后的版本中,minimize 增加了一个location 选项,这个选项有两个选择:true 和false .它的省缺值是false ,如果将location 的值设置为true ,minimize 就可以获得极值点的位置.例如:>minimize((x+1)/(x^2+1),location=true);−√2(−1−√2)2+1,{[{x =−1−√2},−√2(−1−√2)2+1]}242第十章数据处理minimize和maximize命令可以在无约束的条件下求函数的极值,对于有约束条件,可以使用extrema命令来求极值.调用extrema的方法有两种.一种是使用readlib命令,另一种是使用student程序包.extrema命令的一般用法是:extrema(expr,constraints,vars,’s’).其中第二个选项contraints是约束条件,如果没有约束条件,可以用{}来代替.例如:>readlib(extrema):>extrema(a*x^2+b*x+c,{},x);{−14b2−4c aa}约束条件通常是等式,例如:>extrema(x*y*z,x^2+y^2+z^2=1,{x,y,z});{max(−13RootOf(3Z2−1),13RootOf(3Z2−1),0),min(−13RootOf(3Z2−1),13RootOf(3Z2−1),0)}第四个选项是一个未赋值的变量,它可以带回达到极值点的位置. >extrema(x*y,x+y=1,{x,y},’s’);{1 4 }>s;{{x=12,y=12}}§10.4线性优化Maple对于线性优化问题提供了专门的程序包simplex.它包括一些处理线性优化的命令,其中最重要的命令是minimize和maximize.这两个命令的基本用法是:minimize(f,C),其中f 是一个线性表达式,它是目标函数,C是一系列线性约束条件的集合,这些约束条件必须是带有≤或≥的线性不等式,simplex程序包不能处理严格不等式.>with(simplex);Warning,new definition for maximizeWarning,new definition for minimize[basis,convexhull,cterm,define zero,display,dual,feasible,maximize,minimize,pivot, pivoteqn,pivotvar,ratio,setup,standardize]§10.4线性优化243需要注意的是,当我们调用simplex时,它给出了maximize和minimize的新定义.此时maximize和minimize是simplex程序包中的命令,而不再是求函数极值的命令.maximize的用法如下:>maximize(300*x1+400*x2,{2*x1+4*x2<=170,2*x1+2*x2<=150,6*x2<=180});{x1=65,x2=10}maximize还有其它的参数,第三个参数说明变量的类型,可用的值为NONNEGATIVE或UNRESTRICTED.省缺参数为UNRESTRICTED,它的含义是对所有变量都是无约束的.NONNEGATIVE的含义是所有的变量都是非负的.第四个参数必须是一个未赋值的变量,它用来返回对最优化的描述,第五个参数可用来返回变换的情况.>maximize(6*x1+10*x2+8*x3+12*x4,{3*x1+2*x2+2*x4<=60,x1+x2+2*x3+x4<=45,2*x1+4*x4<=90, 4*x1+4*x3+3*x4<=120,x1+2*x2+2*x3+2*x4<=30},NONNEGATIVE,’eqsys’);{x3=0,x2=0,x4=152,x1=15}>eqsys;{x4=−34SL5+152+14SL1−x2−32x3,SL2=452+14SL1+14SL5−32x3,SL3=30+4x2+2SL5+4x3,SL4=752+54SL1+3x2+14SL5−72x3,x1=−12SL1+15+12SL5+x3}需要说明的是maximize在上例中仅得到了一个最优解,还有一个最优解x1=x2=x3=0,x4= 15并没有求出来.。
