概率统计在化工中应用技巧二

 过冷水诚挚邀请你加入Matlab仿真秀官方交流群进行Matlab学习、问题咨询、 Matlab相关资料下载，qq：927550334 

        在采用MATIAB处理各种实验数据时,首先需要将数据输入MATLAB。输入MATLAB的数据可以采用不同的数据类型保存，数值矩阵、单元数组和结构体都可以用于保存实验数据。数值矩阵用于表示实验数据最简单,但只储存一种类型的数据,而且要求矩阵的各维数据数量相等。数据不等时,有时可以使用NaN补足。存储时一般每列元素表示相同条件下的重复数据,每行元素表示实验条件不同时的数据;由于只能表示数值数据,因此很多实验信息则只能采用另外的方法保存。当需要保存实验数据的类型不止一种时,可以采用结构体或单元数组进行储存。这可以将一些相关联的信息存储在一个变量中,方便查找。但是MATIAB大多运算函数仅针对数值元素操作,必须将结构体和单元数组中的数值数据取出才能进行其他运算。总之,各种数据类型储存实验结果各有优势,可以根据具体要求灵活选择。比如如下所示图例

1采用结构体变量表示数据,并计算催化剂组成为A时的平均转化率

Cat(1).comp='A';Cat(1).load='A';Cat(1).ptr='A';Cat(1).R=52.7;
Cat(2).comp='B';Cat(2).load='A';Cat(2).ptr='A';Cat(2).R=60.1;
Cat(3).comp='A';Cat(3).load='B';Cat(3).ptr='A';Cat(3).R=53.5;
Cat(4).comp='B';Cat(4).load='B';Cat(4).ptr='A';Cat(4).R=62.3;
Cat(5).comp='A';Cat(5).load='A';Cat(5).ptr='B';Cat(5).R=55.4;
Cat(6).comp='B';Cat(6).load='A';Cat(6).ptr='B';Cat(6).R=66.0;
Cat(7).comp='A';Cat(7).load='B';Cat(7).ptr='B';Cat(7).R=59.1;
Cat(8).comp='B';Cat(8).load='B';Cat(8).ptr='B';Cat(8).R=68.0;
RT=0;n=0;
for i=1:length(Cat)
    if Cat(i).comp=='A';
        RT=RT Cat(i).R;
        n=n 1;
    end
end
Rava=RT/n
Rava =
   55.1750

2采用数值数组表示实验转化率,单元数组表示制备条件;计算催化剂组成为B时的平均转化率;

Yield=[52.7 60.1 53.5 52.3 55.4 66.0 59.1 68.0];
Composition={'A';'B';'A';'B';'A';'B';'A';'B'};
Load={'A';'A';'B';'B';'A';'A';'B';'B'};
Ptr={'A';'A';'A';'A';'B';'B';'B';'B'};
RT=0;n=0;
for i=1:length(Composition)
    if Composition{i}== 'B'
        RT=RT  Yield(i);
        n=n 1;
    end
end
Rava=RT/n

3采用合适的方法表示实验数据,计算催化剂组成变化对产率的影响,即组成为A时的平均产率与组成为B时的平均产率之差。

ExpData=[1 1 1 52.7;
    -1 1 1 60.1;
    1 -1 1 53.5;
    -1 -1 1 62.3;
    1 1 -1 55.4;
    -1 1 -1 66.0;
    1 -1 -1 59.1;
    -1 -1 -1 68.0];
RD=ExpData(:,1)'*ExpData(:,4)/4

        由以上程序可见,第一种方法采用结构体表示结果赋值比较复杂,但信息表示清楚;第二种表示方法比较简洁,这种方法在MATIAB统计推断问题中经常使用;第三种方法最为简单,而且由于完全是数值矩阵,后续计算问题非常简便,不利之处在于由于对于初始信息进行了一次编码, 增加了犯错误的概率。

样本数字特征的计算

MATLAB计算样本数字特征的函数表格

         mean、geomean、median和 harmmean的使用方法相同,以 mean为例说明。mean函数计算数组不同维上的平均值。 当X为向量时, mean(X)返回X的平均值;如果X为矩阵,则mean将计算X的每列元素的平均值,返回值M为一行向量。

M=mean(X,dim)

此种调用将计算数组指定维 dim 上的平均值，例如 mean(X, 2) 计算X的行向量平均值。

M=Trimmean(X, percent)
M=Trimmean(X, percent, dim)
M=Trimmean(X, percent, fiag)
M=Trimmean(X, percent, fiag, dim)

        输人变量中的X为需要计算的数据,当X是向量时,M为其切尾平均值,当X为二维矩阵﹐则计算每列元素的切尾平均值。percent取值在0~100之间,为需要去掉的样本数据比例,例如当percent为20时, 表示计算时将去掉最大和最小的1括10个元素的向量,则将去掉1个最大值和1个最小值。如果k=length(X)*percent/100/2不等于整数时,可以采用flag 指定的方法取整flag 的取值包括: ' round'表示取最近整数;'floor'取较小整数;' weight'切尾端点有权重的平均值。

     表示样本变异特征的几个函数使用也很简单,最简单的格式是只有一个输入参数,即表示样本数据的向量或矩阵X,函数将计算获得对应的变异特征,例如std(X)表示计算X样本的标准差。其他特殊规定可以通过MATLAB帮助自行学习。

样本数据的图形化表示

        在获得样本数据后,进行图形化表示将有助于辨识数据的质量并进行初步的推断, MATLAB的二维绘图plot命令可以完成图形化工作,除此之外,MATLAB还提供了其他统计可视化函数,以下将介绍散点图、盒状图、误差图和直方图等几种。

散点图

        散点图即将数据点以散点的形式绘制在图形上,这样可以观察数据的大致趋势及有无异常。MATLAB的plot命令可以绘制散点图,不过统计工具箱中的函数gscatter绘制散点图更为有效。gscatter的调用格式如下。

gscatter(x，y，group)
gscatter(x，y，group，clr，sym，siz)

        其中,输入变量α，y表示原始数据,分别对应α轴、y 轴; group为分组信息,可以是向量、字符数组或单元数组,表示了x，y来源的不同; clr、sym和 siz分别表示绘制散点的颜色、符号和大小。

Yield=[52.7 60.1 53.5 52.3 55.4 66.0 59.1 68.0];
Composition={'A';'B';'A';'B';'A';'B';'A';'B'};
gscatter(1:length(Yield),Yield,Composition,'rb','o*')
xlabel('Sample Number')
ylabel('Yield [%]')
box on

Yield=[52.7 60.1 53.5 52.3 55.4 66.0 59.1 68.0];
Load={'A';'A';'B';'B';'A';'A';'B';'B'};
gscatter(1:length(Yield),Yield,Load,'rb','o ')
xlabel('Sample Number')
ylabel('Yield [%]')
box on

盒状图

        样本数据由于存在误差,因此如果重复实验,则数据应随机分布在平均值附近。采用盒状图可以将样本均值及其他分布性质绘制在一起。盒状图的绘制命令为: boxplot

boxpIot(x，group)
boxplot(. . . ，'Name' ,val，. ..)

其中x为需要绘制的数据,可以是矩阵或向量,group表示x的来源的不同,与gscatter中group 的规定相同。'Name'和 val控制盒状图的格式。

试绘制以上数据的盒状图,并根据图形判断两种分析方法是否存在差异。首先绘制图形表示组成的影响,程序如下。

x1=[46.0 46.3 54.2 54.6 50.8 47.3 48.5 36.2];
x2=[51.8 50.0 53.2 49.3 52.8 66.2 52.3 50.5 ];
boxplot([x1',x2'],'labels', {'SPG','AES'},'color','bk')
ylabel('Fe ion concentration [ppm]')

        图中盒子上下边缘和中间直线表示数据的25%分位数(q₁)、75%分位数(q3)和中位数,盒子外的竖线表示了1.5倍q₃-q₁长度,超过这一数值后的数据被称为异常数据,以十字符表示。由上图结果看,两种方法所得结果数据彼此重叠,说明两种方法没有明显区别。具体的统计推断可在定量分析后做出。

误差图

在曲线图上叠加各数据点的误差,则形成误差图。MATLAB的errorbar命令可用于误差图的绘制,其调用格式如下。

errorbar(Y,E)
errorbar(X,Y,E)
errorbar(X,Y,L,U)

        Y表示了y轴数据;E表示误差大小, errorbar将把2*E间距(置信区间)表示在图形上;X为x轴数据;L,U为指定的误差下、上限。在搅拌釜中进行两个配方的对比实验,每个配方进行了三次重复实验,实验获得的反应时间和产品收率如表所示。采用errorbar绘制下表数据,置信区间间距为3倍的样本标准差。根据图形能否确定两个配方的优劣?

Pro.A.tim=[1.1;2.2 ;3.2;4.1]
Pro.A.test=[10.1,15.4,22.5,31.6;11.5,12.6,21.6,32.9;9.6,17.8,22.9,32.4;]'
Pro.B.tim=[0.8;1.5;2.9;4.0];
Pro.B.test=[8.5,14.0,21.4,35.1;9.6,15.5,22.6,33.0;8.8,15.2,22.0,31.0]'
Amean=mean(Pro.A.test');Bmean=mean(Pro.B.test');
AError=std(Pro.A.test');BError=std(Pro.A.test');
figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
errorbar(Pro.A.tim',Amean,AError,'MarkerFaceColor',[1 0 0],'Marker','.','LineStyle',':','LineWidth',3)
hold on
errorbar(Pro.B.tim',Bmean,BError,'MarkerFaceColor',[0 1 0],'Marker','.','LineWidth',3)
ylabel('Yield [%]','FontName','Times New Roman');
xlabel('Reaction Time [min]','FontName','Times New Roman');
legend('Composition A','Composition B')
set(axes1,'FontName','Times New Roman','FontSize',15,'LineWidth',3);

绘制结果如图所示。由此可见,虽然方法B的最终收率看似略高于A,但由于置信区间包括了配方A的数据,因此不能可靠地判断B配方比A配方更好。

        过冷水发表于 仿真秀 平台原创文章，未经授权禁止私自转载，如需转载请需要和作者沟通表明授权声明，未授权文章皆视为侵权行为，必将追责。如果您希望加入Matlab仿真秀官方交流群进行Matlab学习、问题咨询、 Matlab相关资料下载均可加群：927550334。

精品回顾

 matlab绘制农夫过河动态图

分子动力学的原子空间运动轨迹演示编程

过冷水带你用matlab制作演示动画

python批量移动文件&重命名代码分享

过冷水和你分享 matlab读取存储各种文件的方法 文末有独家金曲分享