matlab用于试验设计回归分析实验结果的例子

本来是转发代码的，后来发现问类似问题的人太多了，干脆根据自己的感想，总结成一篇原创吧。

快速入门的方法，软件推荐

希望不要再来问我了。——建议自己先去学一款适用于DOE的软件，同学们！（我自己也是主要靠软件、而不是靠学习统计原理来做事的）

这些软件如下：

1. Minitab, 大部分搞six sigma, lean，DMAIC/DFSS的公司采购的都是这个软件；(流行，但不见得强大，公司用户的统计知识水平可想而知）

2. Statistica , 比Minitab更强大、使用方法也相似（不知道跟minitab之间是谁学了谁），以前说四大数学软件，我一直误以为是matlab,mathematica,maple和statistica

3. SAS JMP 这SAS应该大名鼎鼎了；该公司还开发了专门搞试验设计的JMP, 以前使用过，感觉：强大，DOE的话，估计比minitab强；

4. Design Expert, 我没有用过，但是，下面一本英文书里推荐的三款软件之一就是它(JMP, Minitab),小巧强大；

5. unscramber也是我在网络上看过的一款类似软件，没有用过，说实话，一款软件强大的足矣了，用不着每种都学：想起那个典故： 上衣口袋别一支钢笔，是中学生，两只是大学生，别三支以上钢笔……是修钢笔的。

这些软件大多有文档和帮助，网络上也有案例，建议全程学习一两个案例之后，搞懂了再深入琢磨。有一个全局的掌控之后才明白哪里该侧重用力。此外

做实验设计之前，最好先搞定这些准备工作：

1. 明白自己实验的目的，只是找显著因子、还是要摸索最佳实验条件？

很多情况下，严格意义的复杂的试验设计可能并不是必要的，有时候可能只需要一些简单的对比、少量序贯试验，或者边试验边调整某些因子就能更轻松达到目的的，强烈不建议还要使用复杂的试验设计。

此外，是不是用复杂的试验设计，还取决于，试验的体系对于自己来说是完全的黑箱、还是灰箱、甚至白箱。也就是，完全两眼一抹黑，毫不了解(黑箱，深度灰箱)，DOE的帮助会是较大的。有些体系，尤其是某些已经被前人发表的论文研究得非常透彻的化学反应体系，其最佳配方、投料比和最优实验条件之类的问题，往往是很清楚的，只须简单验证即可，DOE的意义就不大了。——所以，开始做之前，尝试通过查阅足够多的文献了解自己的研究对象，尽可能获取详细的信息，其实是DOE之前必须的、很重要的准备工作。有时候可能找相关的专家、前辈的老师或师兄师姐给点参考意见，往往也比直接做DOE更能起到事半功倍效果。所以，这个不可以不知道。

不过我知道，很多情况下，做一个复杂点的试验设计，不仅仅为了某个最佳的条件，还为了一个漂亮的分析结果、方便发论文、混毕业……这种情况下，可以理解。

2. 清楚自己实验的成本、时间、人力等局限性，以及自己可以接受的投入跟目的之间的是否有了恰当的平衡？

在某某公司的时候，某种试验的成本相当昂贵，设备都是几千万美元的设备，人也都是昂贵的人，材料都是不可回收的昂贵材料，但Six sigma是一种类似政治口号的东西，大家谁用得不漂亮就不受尊重，升迁和加薪的机会就渺茫，所以，即使成本昂贵，也在被反复使用，然后得到一些未必有太大价值的结果，或者不用这样的方案也能用廉价方式得到的结果。所以，有些时候，不得不承认，老板被six sigma的咨询顾问公司忽悠之后(这老板在世界顶级的公司都待过，造过飞机发动机和飞机），six sigma运动化，半懂不懂地强推统计和数据处理方法，对一个公司的影响是绝对负面的。——好在百足之虫，经得起这样的折腾；好在听说后来的大老板是个工科出生的，一上来就把这种运动的影响逐渐消除掉、淡化了。

3.  根据实验的目的、约束条件、确定试验设计方案和分析方法。

全析因试验设计通常是笨拙的、穷尽各种可能，形成一个全排列，实验量很大。但是在有些场合，我自己经历过的，还是可行的，因为每个试验点成本很低、可以轻易并行地大规模同步实验，即使每个点都有重复、多至上千个试验点，工作量也在当时的实验目的及预算所能轻松承担的范围之内（所以我选的是带有重复的全析因试验设计方案）；而且，只需要作一个方差分析，找显著因子就可以大功告成了，没有找最优实验条件之类的要求（所以，只须用Minitab中的方差分析就可以了；当然更加可以用回归分析，而且几乎不会增加什么特别的成本，但是，已经不需要了；）。

——但一般情况下，实验成本高、周期长的时候，才考虑其它方案。国外用response surface （回归分析找最优点）的比较多；中日东亚的确对Taguchi ,田口、正交试验设计更加情有独钟；此外，还有实验成本更低的均匀实验设计（因为试验次数更少，通常回归分析效果更好）。

4. 根据前面的理解，选择合适的工具软件。

以往那种查表格、手工计算的方式效率低、容易出错，合适的软件： Minitab、SAS JMP、unscramber、Design Expert 等 都是不错的选项；其它统计软件也通常都能完成类似这样或那样的任务；——即使前面特别提到的这些软件，以及其它，并不是所有的软件都提供所有这些已经发表了的方法的表格生成和结果分析的功能。所有这些软件也都是有自己的侧重和局限性的。使用的时候要先弄清自己需要的是什么，然后选合适的工具。

5. 当常规的试验设计概念或软件不能满足条件时候怎么办？ 下面的这个是，当一般的试验设计概念无法满足要求的情况下，如何计算实验表格以及分析结果。

详情可以参考 方开泰、马长兴， 科学出版社 ：正交与均匀试验设计 一书的理论解释。

其实推荐用于D-optimal试验设计表格生成的函数和代码，不懂的地方请参考matlab相关函数的对应的文档。其它大型计算类软件应该也能完成类似的功能，不再提。

为什么突然选择了matlab? 因为有现成的好用的函数。

matlab对正交试验设计的支持是很弱的，只有简单的全析因和部分析因、响应曲面法和D-最优设计几种可选

想要实现正交表格，还须自力更生，利用正交表实际是特殊的D-最优表格的事实。（参考 方开泰、马长兴：正交和均匀试验设计）

代码实现在这里

但是一点也不推荐用matlab来做正交试验设计表格，因为这不是它的专长。

https://icme.hpc.msstate.edu/mediawiki/index.php/DOE_with_MATLAB_3

DOE with MATLAB 3

Jump to: navigation, search

Abstract

This example shows how to do multilevel full factorial designs and Taguchi designs using MATLAB.

Author(s): Mark A. Tschopp

Introduction to D-optimal Designs in MATLAB

From MATLAB help:

"Traditional experimental designs (Full Factorial Designs, Fractional Factorial Designs, and Response Surface Designs) are appropriate for calibrating linear models in experimental settings where factors are relatively unconstrained in the region of interest. In some cases, however, models are necessarily nonlinear. In other cases, certain treatments (combinations of factor levels) may be expensive or infeasible to measure. D-optimal designs are model-specific designs that address these limitations of traditional designs.

A D-optimal design is generated by an iterative search algorithm and seeks to minimize the covariance of the parameter estimates for a specified model. This is equivalent to maximizing the determinant D = |XTX|, where X is the design matrix of model terms (the columns) evaluated at specific treatments in the design space (the rows). Unlike traditional designs, D-optimal designs do not require orthogonal design matrices, and as a result, parameter estimates may be correlated. Parameter estimates may also be locally, but not globally, D-optimal."

Multilevel Full Factorial Design

What if we have four factors (a, b, c, d) at three levels (1 2 3)?

nfactors = 4; nlevels = 3; nruns = nlevels^nfactors; % 81 runs dfF = sortrows(rowexch(nfactors,nruns,'l','cat',1:nfactors,'levels',nlevels*ones([1,nlevels]),'tries',1))

which gives...

dfF = 1 1 1 2 1 1 1 2 1 1 2 1 1 1 2 2 1 1 2 3 1 1 2 3 1 1 2 3 1 1 3 1 1 1 3 3 1 2 1 1 1 2 1 3 1 2 1 3 1 2 2 1 1 2 2 2 1 2 3 1 1 2 3 1 1 2 3 1 1 2 3 2 1 3 1 1 1 3 1 2 1 3 1 3 1 3 1 3 1 3 2 1 1 3 2 3 1 3 3 2 1 3 3 2 1 3 3 2 2 1 1 2 2 1 1 2 2 1 1 3 2 1 2 1 2 1 2 2 2 1 2 3 2 1 3 1 2 1 3 2 2 1 3 2 2 2 1 1 2 2 1 1 2 2 1 2 2 2 1 2 2 2 2 2 2 2 2 3 2 2 3 3 2 2 3 3 2 2 3 3 2 3 1 1 2 3 1 3 2 3 2 1 2 3 2 1 2 3 2 1 2 3 2 2 2 3 3 1 2 3 3 3 2 3 3 3 3 1 1 1 3 1 1 1 3 1 1 1 3 1 1 3 3 1 2 1 3 1 3 1 3 1 3 2 3 1 3 3 3 1 3 3 3 2 1 2 3 2 1 2 3 2 2 1 3 2 2 2 3 2 2 2 3 2 2 3 3 2 2 3 3 2 3 1 3 2 3 3 3 3 1 2 3 3 1 3 3 3 1 3 3 3 2 1 3 3 2 2 3 3 2 3 3 3 3 1 3 3 3 2 3 3 3 2

What if we have four factors (a, b, c, d) at mixed levels (1 2 for 'a', 1 2 3 for 'b','c','d')? Just use...

nfactors = 4; nruns = 54; % 2*3*3*3 runs dfF = sortrows(rowexch(nfactors,nruns,'l','cat',1:nfactors,'levels',[2 3 3 3],'tries',1))

Taguchi design

Here are the commands for L8, L12, L25 Taguchi design arrays. Notice that you can used the 'bounds' command to change the factor levels from 1 and 2 to evenly spaced numbers between bounds (-1 and +1 in this case). However, for the five factor levels in the L25 design, I removed the command to show the factor levels as 1-5. Within the Taguchi designs, 7 factors at 2 levels can be used in the L8 design, 11 factors at 2 levels can be used in the L12 design, and 5 factors at 5 levels can be used in the L25 design.

i = 7; L8 = sortrows(rowexch(i,8,'l','cat',1:i,'bounds',[-1*ones([1,i]);ones([1,i])],'levels',2*ones([1,i]),'tries',100)); i = 11; L12 = sortrows(rowexch(i,12,'l','cat',1:i,'bounds',[-1*ones([1,i]);ones([1,i])],'levels',2*ones([1,i]),'tries',100)); i = 5; L25 = sortrows(rowexch(i,25,'l','cat',1:i,'levels',5*ones([1,i]),'tries',100));

which give the following arrays...

L8 array

Again, seven factors at two levels (-1 and +1).

L8 = -1 -1 -1 1 1 1 -1 -1 -1 1 -1 1 -1 1 -1 1 -1 -1 -1 1 1 -1 1 1 1 -1 -1 -1 1 -1 -1 -1 -1 -1 -1 1 -1 1 1 -1 1 1 1 1 -1 1 1 -1 1 1 1 1 -1 1 1 -1

L12 array

Again, eleven factors at two levels (-1 and +1).

L12 = -1 -1 -1 1 -1 -1 1 1 -1 -1 1 -1 -1 -1 1 1 1 -1 -1 1 -1 -1 -1 -1 1 -1 1 -1 -1 1 -1 1 -1 -1 1 -1 -1 -1 1 1 1 1 1 -1 -1 1 1 -1 1 -1 1 -1 1 -1 1 -1 1 1 1 -1 1 -1 -1 -1 1 1 1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 -1 1 -1 -1 1 1 -1 -1 -1 -1 1 -1 1 1 1 1 1 1 1 1 1 1 1 -1 -1 1 1 -1 1 -1 -1 1 1 1 -1 1 1 -1 1 -1 -1 1 -1 1 1 1 1 -1 -1 -1 1 1 -1 -1

L25 array

Again, five factors at five levels (1-5).

L25 = 1 1 2 5 3 1 2 4 1 2 1 3 5 2 4 1 4 3 4 5 1 5 1 3 1 2 1 4 4 4 2 2 1 5 5 2 3 3 3 2 2 4 2 2 1 2 5 5 1 3 3 1 3 1 1 3 2 2 3 4 3 3 1 4 3 3 4 5 5 2 3 5 4 2 5 4 1 5 3 5 4 2 3 2 3 4 3 4 5 1 4 4 1 1 4 4 5 2 4 2 5 1 1 2 2 5 2 5 4 1 5 3 2 1 5 5 4 4 3 3 5 5 3 5 4

General

You can use these functions (sortrows, rowexch) to create other d-optimal designs. In some cases, obtaining the d-optimal designs for Taguchi matrices may not be very efficient with these functions or may require too much memory. In these cases, you can always enter the design matrix directly from Taguchi's generated matrices.

In general, if you want a different number of factors and factor levels than the traditional Taguchi design arrays, you can produce them using the following techniques. It is recommended that your designs are balanced. That is, the number of levels for each factor (each column above) are equal. If you have 8 runs with two factor levels for factor 'a', it is best to have a design where there are four combinations with 'a' at a low level and four combinations with 'a' at a high level, etc.

Go Back

对3因素混合水平[3,3,4]和有交互作用,采用18次试验的方案,D-optimal的试验设计表格如下:

warning off;
sortrows(rowexch(3,18,'interaction','categorical',1:3,'levels',[3,3,4],'tries',1000,'maxiter',50))

方差分析用matlab的anovan就可以了;

回归分析的一个例子matlab代码来自华东理工大学黄华江博士的书

function DOE 
% D-优化试验设计 
% 
%   Author: HUANG Huajiang 
%   Copyright 2003 UNILAB Research Center,  
%   East China University of Science and Technology, Shanghai, PRC 
%   $Revision: 1.0 $  $Date: 2003/07/12 $ 
 
clear all 
clc 
 
% 各因素的水平 
p = [470, 300,120 
     285, 190, 65 
     100, 80, 10]; 
 
% D-优化试验设计：先生成settings，再将settings转变为相应实验条件expCond 
settings = cordexch(3,13,'q');      % settings: 因子设置矩阵 
mr = p(2,:);                        % mr: middle row, i.e., middle level 
mr = mr(ones(13,1),:); 
hr = (p(1,:) - p(3,:))/2; 
hr = hr(ones(13,1),:); 
expCond = settings.*hr + mr;        % expCond: settings的相应实验条件 
 
% 由反应模拟器生成实验数据data 
p1 = expCond(:,1); 
p2 = expCond(:,2); 
p3 = expCond(:,3); 
y = zeros(13,1); 
for k = 1:13 
    y(k)  = 1.25*(p2(k) - p3(k)/1.5183)./(1 + 0.064*p1(k)     ... 
            + 0.0378*p2(k) + 0.1326*p3(k))*normrnd(1,0.02); 
end 
data = [expCond y]     % data为实验数据矩阵（实验条件及其对应的反应速率） 
 
% 数据分析: 由非线性模型Nonlinear Model估计参数 
x = data(:,1:3);  
y = data(:,4); 
xname = str2mat('Hydrogen','n-Pentane','Isopentane'); 
yname = 'Reaction Rate'; 
beta0 = [1.2 0.1 0.01 0.1 1.5];             % 参数初值 
nlintool(x,y,@hougen,beta0,[],xname,yname); 
% rstool(x,y,[],[],xname,yname); 
rstool(x,y,'quadratic',[],xname,yname); 
 
[beta,resid,j] = nlinfit(x,y,@hougen,beta0) 
ci = nlparci(beta,resid,j) 
 
% 参数辨识结果:β1、β2、... β5 
fprintf('Estimated Parameters:\n') 
fprintf('\tβ1 = %.2f ± %.2f\n',beta(1),ci(1,2)-beta(1)) 
fprintf('\tβ2 = %.3f ± %.3f\n',beta(2),ci(2,2)-beta(2)) 
fprintf('\tβ3 = %.4f ± %.4f\n',beta(3),ci(3,2)-beta(3)) 
fprintf('\tβ4 = %.4f ± %.4f\n',beta(4),ci(4,2)-beta(4)) 
fprintf('\tβ5 = %.4f ± %.4f\n',beta(5),ci(5,2)-beta(5)) 
 

发现还有同学有下面这样的问题，特地转过来：

为何 F 值和 p 值估计值在输出中显示为星号？

星号表示因模型饱和且误差没有足够的自由度而无法计算的缺失值。

假设一个饱和全因子 DOE 模型的示例：包含因子 A、B 和 C 而无仿行、无中心点且无区组的 3 因子、两水平设计。此设计将进行 8 次试验。

分析设计时，选择通过包含所有主效应（A、B、C）和所有交互作用项（AB、AC、BC、ABC）来拟合饱和模型。生成的方差分析表会用星号表示残差误差的 SS 值、残差误差的 MS 值、所有 F 统计量和所有 p 值：

方差分析   

来源    自由度    AdjSS    AdjMS    F值    P值
模型    7    71.9880    10.2840    *    *    
线性    3    38.9547    12.9849    *    *    
C5    1    7.0882    7.0882    *    *    
C6    1    7.9818    7.9818    *    *    
C7    1    23.8848    23.8848    *    *    
2因子交互作用    3    32.7537    10.9179    *    *    
C5*C6    1    12.0209    12.0209    *    *    
C5*C7    1    6.5509    6.5509    *    *    
C6*C7    1    14.1818    14.1818    *    *    
3因子交互作用    1    0.2796    0.2796    *    *    
C5*C6*C7    1    0.2796    0.2796    *    *    
误差    0    *    *    
合计    7    71.9880

缺失值在表格中是因为 Minitab 不可能计算这些统计量。由于残差误差有 0 个自由度 (DF)，因此不可能计算这些统计量，如以下计算所示：

• 总自由度 = 运行次数 - 1
• 主效应自由度 = 因子水平数 - 1
• 交互作用效应自由度 = 分量因子的自由度相乘
• 残差误差自由度 = 总自由度 - 模型中包含的所有项的自由度之和

因此，使用先前的示例：

• 总自由度 = 8 - 1 = 7（8 行数据）
• 因子 A 的自由度 = 2 - 1 = 1（因子 A 有 2 个水平）
• 因子 B 的自由度 = 2 - 1 = 1
• 因子 C 的自由度 = 2 - 1 = 1
• 交互作用 AB 的自由度 = (1)*(1) = 1（因子 A 有 1 个自由度，因子 B 有 1 个自由度）
• 交互作用 AC 的自由度 = (1)*(1) = 1
• 交互作用 BC 的自由度 = (1)*(1) = 1
• 交互作用 ABC 的自由度 = (1)*(1)*(1) = 1
• 残差误差的自由度 = 7 - (1 + 1 + 1 + 1 + 1 + 1 + 1) = 0

误差的自由度为零会导致计算失败，如下所示。通过将 Adj SS 列中的值除以 DF 列中的对应值来计算 Adj MS 列中的每个值（因子 A 的 Adj MS = Adj SS/DF = 0.0621/1 = 0.0621）。但是，由于不可能将任何值除以 0 个自由度，因此，无法计算残差误差的 Adj MS（通常称为均方误 (MSE)）。

而且，Minitab 通过将每个 Adj MS 值除以 MSE 来计算表格 F 列中的每个值。例如，因子 A 的 F 值等于 0.0621/MSE。但是由于无法计算 MSE，因此也无法计算 F 值。

最后，根据 F 统计量计算 p 值。因此，如果 F 缺失，则 p 值也必须缺失。

当存在包含一个仿行的两水平设计且模型中包含所有项时，方差分析表中会出现缺失的 p 值和 F 统计量。要补救这种情况，请重新拟合不含一个或多个交互作用项的模型。为确定要从饱和模型中删除的最高阶交互作用，请使用效应图估计交互作用的统计显著性。

例如，如果选择统计 > DOE > 因子 > 分析因子设计，单击模型按钮，并从模型中删除 ABC 交互作用项，则 Minitab 可以为主效应和双因子交互作用计算方差分析表中的所有值。

方差分析来源 自由度 Adj SS Adj MS F 值 P 值模型 6 71.7084 11.9514 42.75 0.117 线性 3 38.9547 12.9849 46.45 0.107 C5 1 7.0882 7.0882 25.35 0.125 C6 1 7.9818 7.9818 28.55 0.118 C7 1 23.8848 23.8848 85.44 0.069 2 因子交互作用 3 32.7537 10.9179 39.05 0.117 C5*C6 1 12.0209 12.0209 43.00 0.096 C5*C7 1 6.5509 6.5509 23.43 0.130 C6*C7 1 14.1818 14.1818 50.73 0.089误差 1 0.2796 0.2796合计 7 71.9880

现在，由于误差还剩余 1 个自由度（这意味着 Minitab 可以计算 MSE、F 和 p 值），因此，Minitab 会计算所有值。