Showing posts with label Matlab. Show all posts
Showing posts with label Matlab. Show all posts

Monday, May 4, 2015

Convert MATLAB history.m to History.xml

Link.


fid = fopen('history.m');

% Create the document node and the root element, 'history'
docNode = com.mathworks.xml.XMLUtils.createDocument('history');
% Get the history node
history = docNode.getDocumentElement;
% Read the first line of the old history.m file
line = fgetl(fid);

% Loop through every line of the old history.m file while there are lines
while ischar(line)
% Add a command element for each command
command = docNode.createElement('command');
% Add the command itself as the child text node
command.appendChild(docNode.createTextNode(line));
% Add the command node back as a child of the command node
history.appendChild(command);
% Get the next line of the old history.m file
line = fgetl(fid);
end

% Close the file
fclose(fid);

% Write the XML file
xmlwrite('History.xml',docNode);
Read more ...

Wednesday, January 14, 2015

Overlay Line Plot on Bar Graph Using Different Y-Axes and change the properties of Bar and Line

For Overlay Line Plot on Bar Graph Using Different Y-Axes. See here.

days = 0:5:35;
conc = [515,420,370,250,135,120,60,20];
temp = [29,23,27,25,20,23,23,27];

[AX,hBar,hLine] = plotyy(days,temp,days,conc,'bar','plot'); %u can switch order of 'bar' and 'plot'.

%Change the Axis properties
set(AX(1),'ycolor','r')
set(AX(2),'ycolor','b')

%Change Line properties
set(hLine, 'color', r)

%Change Bar properties
set(hBar, 'facecolor', r)
set(hBar, 'edgecolor', g)

%Transparent bar
ch = get(hBar, 'child')
set(ch, 'facea', 0.5)
Read more ...

Monday, September 15, 2014

Sunday, September 15, 2013

Java 1.6.0_51 breaks MATLAB 2012b

The was an issue in the Java security updates that Apple released for Mac OS X (as mentioned above). These updates install Java 1.6.0_51 and the exact build number looks like this (ending with M4508):
Java™ SE Runtime Environment (build 1.6.0_51-b11-456-10M4508) You can check if you have installed the xM4508 versions by running any one of the following commands on the Terminal of your Mac OS X:
java -version 
OR
/usr/libexec/java_home -v 1.6 -exec java -version 
If you have previously installed the xM4508 versions of the Java updates you can upgrade to the fixed xM4509 version by manually installing the following update:
Mac OS 10.7.x and 10.8.x:
================
Mac OS 10.6.x:
================
Confirm that you have the updated Java version by executing any of the above two commands. The exact build number looks like this (ending with M4509):
Java(TM) SE Runtime Environment (build 1.6.0_51-b11-457-10M4509)
updates already (ending with M4508), you should receive this fix (ending with M4509) automatically through Software Update / App Store.

source: here
Read more ...

Friday, June 21, 2013

MATLAB FUNCTION accumarray

A = accumarray(subs,val,sz,fun,fillval)
  sub:提供累计信息的指示向量
   val:提供累计数值的向量
    sz:控制输出向量A的size
   fun:用于计算累计后向量的函数,默认为@sum,即累加
fillval:填补A中的空缺项,默认为0


例子:

  1. myfun = @(x) 1./(2*numel(x)) * sum(x);  

myfun = 

    @(x)1./(2*numel(x))*sum(x)

  1. subs=[7,2,5,7,2]';  

subs =

     7
     2
     5
     7
     2

  1. val=1:5;  

val =

     1     2     3     4     5

  1. sz=[10,1];  

sz =

    10     1

  1. A=accumarray(subs,val,sz,myfun,nan);  

A =

       NaN
    1.7500
       NaN
       NaN
    1.5000
       NaN
    1.2500
       NaN
       NaN
       NaN

个人对此函数的理解:
首先,函数根据参数中sz,生成一个中间矩阵B,此例中sz=[10,1],所以B是一个10*1的矩阵。
然后,函数根据subs中的指示,将val中的数值摆放到B中。
  subs  val
     7     1
     2     2
     5     3
     7     4
     2     5
意思就是,val中的1被扔到了B中第7的位置,同样,2扔到B中第2个位置,3到B中5,4到B中7,5到2。
这样B中的2,5,7位就被摆放上了数值,而且,2和7两个位置上被摆上了2个数值。

B =

   NaN
   2,5
   NaN
   NaN
   3
   NaN
   1,4
   NaN
   NaN
   NaN
其他还有7个位置并没有被摆放任何数值,我们暂时不管他,给他扔上NaN。
然后,accumarray函数就会根据自定义的函数myfun,来计算了。拿什么算呢,拿的就是这个中间矩阵B,我们定义的myfun意义就是均值的一般,所以也可以定义的时候写成:

  1. myfun = @(x) mean(x)./2;  

这下就简单了,就相当于mean(B)./2,然后得出的结果就是A了。
A =

       NaN
    1.7500
       NaN
       NaN
    1.5000
       NaN
    1.2500
       NaN
       NaN
       NaN
个人理解大概是个这样的过程,接下来再去看matlab的help文件就比较容易懂了。
Read more ...

MATLAB中的共轭转置与非共轭转置

对已知矩阵A,MATLAB为我们提供了两种转置运算。
A.' 非共轭转置

A' 共轭转置

单纯地共轭用:conj()
非共轭可以用:transpose()

example:
  a =
        12.0000                  0 + 2.0000i         5.0000         
        0                             5.0000               4.0000 
>> a'
             ans =
                      12.0000                  0         
                      0 - 2.0000i              5.0000         
                      5.0000                    4.0000         
>> a.'
           ans =

                   12.0000                  0         
                  0 + 2.0000i              5.0000         
                  5.0000                    4.0000

Read more ...

MATLAB 中用 spline 代替 solve/finvers 数字求解反函数

example:

  1. Cg=sym('exp(-((-log(u))^alpha + (-log(v))^alpha)^(1/alpha))');  
  2.   
  3. cu=diff(Cg,'u');  
  4. c=subs(cu,{'u','alpha'},[0.8,2]);  
  5.   
  6. vr=linspace(0,1,1000);  
  7. tr=subs(c,'v',vr);  
  8.   
  9. t=rand(100,1);  
  10. v=spline(tr,vr,t);  

如果直接用solve/finvers求解c。matlab提示无解。
##
cu是Cg对u的偏导,c为偏导函数中带入u=0.8 alpha=2得到的只含变量v的函数。
vr和tr是为了生成一定数量的pairs,以便spline函数求反函数的时候插值。

  1. plot(vr,tr,'--')  

可以看到c函数的图。
t为随机生成的100个点,v是这100个点对应的反函数值!
如果是:

  1. v=spline(vr,tr,t);  

那么就是简单的插值。
help spline
可以看到spline的相关功能和用法。
Read more ...

利用 MATLAB 对文件批量重命名

>> matlab 中 strrep 函数可以更改文件扩展名。
>> matlab 中 unix system ! 都可以让你在 matlab 中执行系统命令。
比如:

  1. !mv a.tex b.tex  
  2. unix('mv a.tex b.tex')  
  3. system('mv a.tex b.tex')  

>> 再加上 dir 函数 以及一个 for-loop 就可以对文件进行批量重命名了。
当然,windows 下应该是用 rename 重命名。
Read more ...

MATLAB 常用的基本数学函数

一、MATLAB常用的基本数学函数
abs(x):纯量的绝对值或向量的长度
angle(z):复数z的相角(Phase angle)
sqrt(x):开平方
real(z):复数z的实部
imag(z):复数z的虚部
conj(z):复数z的共轭复数
round(x):四舍五入至最近整数
fix(x):无论正负,舍去小数至最近整数
floor(x):地板函数,即舍去正小数至最近整数
ceil(x):天花板函数,即加入正小数至最近整数
rat(x):将实数x化为分数表示
rats(x):将实数x化为多项分数展开
sign(x):符号函数 (Signum function)。
当x<0时,sign(x)=-1
当x=0时,sign(x)=0
当x>0时,sign(x)=1
rem(x,y):求x除以y的馀数
gcd(x,y):整数x和y的最大公因数
lcm(x,y):整数x和y的最小公倍数
exp(x):自然指数
pow2(x):2的指数
log(x):以e为底的对数,即自然对数或
log2(x):以2为底的对数
log10(x):以10为底的对数

二、MATLAB常用的三角函数
sin(x):正弦函数
cos(x):馀弦函数
tan(x):正切函数
asin(x):反正弦函数
acos(x):反馀弦函数
atan(x):反正切函数
atan2(x,y):四象限的反正切函数
sinh(x):超越正弦函数
cosh(x):超越馀弦函数
tanh(x):超越正切函数
asinh(x):反超越正弦函数
acosh(x):反超越馀弦函数
atanh(x):反超越正切函数

三、适用於向量的常用函数
min(x): 向量x的元素的最小值
max(x): 向量x的元素的最大值
mean(x): 向量x的元素的平均值
median(x): 向量x的元素的中位数
std(x): 向量x的元素的标准差
diff(x): 向量x的相邻元素的差
sort(x): 对向量x的元素进行排序(Sorting)
length(x): 向量x的元素个数
norm(x): 向量x的欧氏(Euclidean)长度
sum(x): 向量x的元素总和
prod(x): 向量x的元素总乘积
cumsum(x): 向量x的累计元素总和
cumprod(x): 向量x的累计元素总乘积
dot(x, y): 向量x和y的内积
cross(x, y): 向量x和y的外积

四、MATLAB的永久常数
i或j:基本虚数单位
eps:系统的浮点(Floating-point)精确度
inf:无限大, 例如1/0
nan或NaN:非数值(Not a number),例如0/0
pi:圆周率 p(= 3.1415926...)
realmax:系统所能表示的最大数值
realmin:系统所能表示的最小数值
nargin: 函数的输入引数个数
nargin: 函数的输出引数个数

五、MATLAB基本绘图函数
plot: x轴和y轴均为线性刻度(Linear scale)
loglog: x轴和y轴均为对数刻度(Logarithmic scale)
semilogx: x轴为对数刻度,y轴为线性刻度
semilogy: x轴为线性刻度,y轴为对数刻度

六、plot绘图函数的叁数
字元 颜色 字元 图线型态
y 黄色 . 点
k 黑色 o 圆
w 白色 x x
b 蓝色 + +
g 绿色 * *
r 红色 - 实线
c 亮青色 : 点线
m 锰紫色 -. 点虚线
-- 虚线

七、注解
xlabel('Input Value'); % x轴注解
ylabel('Function Value'); % y轴注解
title('Two Trigonometric Functions'); % 图形标题
legend('y = sin(x)','y = cos(x)'); % 图形注解
grid on; % 显示格线

八、二维绘图函数
bar 长条图
errorbar 图形加上误差范围
fplot 较精确的函数图形
polar 极座标图
hist 累计图
rose 极座标累计图
stairs 阶梯图
stem 针状图
fill 实心图
feather 羽毛图
compass 罗盘图
quiver 向量场图
Read more ...

MATLAB FUNCTION fmincon

x = fmincon(fun,x0,A,b)
x = fmincon(fun,x0,A,b,Aeq,beq)
x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub)
x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)
x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options,P1,P2, ...)
[x,fval] = fmincon(...)
[x,fval,exitflag] = fmincon(...)
[x,fval,exitflag,output] = fmincon(...)
其中,x, b, beq, lb,和ub为线性不等式约束的上、下界向量, A 和 Aeq 为线性不等式约束和等式约束的系数矩阵矩阵,fun为目标函数,nonlcon为非线性约束函数。
显然,其调用语法中有很多和无约束函数fminunc的格式是一样的,其意义也相同,在此不在重复介绍。对应上述调用格式的解释如下:
x = fmincon(fun,x0,A,b) 给定初值x0,求解fun函数的最小值x。fun函数的约束条件为A*x <= b,x0可以是标量或向量。
x = fmincon(fun,x0,A,b,Aeq,beq) 最小化fun函数,约束条件为Aeq*x = beq 和 A*x <= b。若没有不等式线性约束存在,则设置A=[]、b=[]。
x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub) 定义设计变量x的线性不等式约束下界lb和上界ub,使得总是有lb <= x <= ub。若无等式线性约束存在,则令Aeq=[]、beq=[]。
x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon) 在上面的基础上,在nonlcon参数中提供非线性不等式c(x)或等式ceq(x)。 fmincon函数要求c(x) <= 0且ceq(x) = 0。
x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options) 用options参数指定的参数进行最小化。
x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options,P1,P2,...) 将问题参数P1, P2等直接传递给函数fun和nonlin。若不需要这些变量,则传递空矩阵到A, b, Aeq, beq, lb, ub, nonlcon和 options。
[x,fval] = fmincon(...) 返回解x处的目标函数值到fval。
[x,fval,exitflag] = fmincon(...) 返回exitflag参数,描述函数计算的有效性,意义同无约束调用。
[x,fval,exitflag,output] = fmincon(...) 返回包含优化信息的输出参数output。
非线性不等式约束nonlcon的定义方法
该参数计算非线性不等式约束c(x)<=0 和非线性等式约束ceq(x)=0。 nonlcon 参数是一个包含函数名的字符串。该函数可以是M文件、内部文件或MEX文件。它要求输入一个向量x,返回两个变量—解x处的非线性不等式向量c和非线性等式向量ceq。例如,若nonlcon='mycon',则M文件mycon.m须具有下面的形式:
function [c,ceq] = mycon(x)
c = ...     % 计算x处的非线性不等式。
ceq = ...   % 计算x处的非线性等式。
若还计算了约束的梯度,即options = optimset('GradConstr','on')
则nonlcon函数必须在第三个和第四个输出变量中返回c(x)的梯度GC和ceq(x)的梯度Gceq。
function [c,ceq,GC,GCeq] = mycon(x)
 c = ...          % 解x处的非线性不等式。
 ceq = ...        % 解x处的非线性等式。
 if nargout > 2   % 被调用的nonlcon函数,要求有4个输出变量。
    GC = ...      % 不等式的梯度。
    GCeq = ...    % 等式的梯度。
 end
4.1应用举例
已知某设计问题可以简化为如下数学模型:
 
显然,此模型属于一个二维约束优化问题。应用fmincon函数求解此优化模型,需要如下几个步骤:
1)编制目标函数的M文件
 在Matlab主窗体的命令行中键入:“edit myobj.m”,并在打开的窗口中编制代码创建目标函数M文件:
function f=myobj(x)
f=2*x(1)^2+2*x(2)^2-2*x(1)*x(2)-4*x(1)-6*x(2);
将其保存为myobj.m备用。
2)编制非线性约数函数的M文件
 若有非线性约束,则应用如下步骤创建约束函数M文件:在Matlab主窗体的命令行中键入:“edit mycon.m”并在打开的窗口中编制相应的代码创建约束函数M文件:
function [c,ceq]=mycon(x)
% 非线性不等式约束条件的表达式,c(1)=...,c(2)=...
c(1)=x(1)+5*x(2)^2-5;
%非线性等式约束条件的表达式
ceq=[];
本例中没有非线性约束,故可以用上述表达方式,也可省略这一步。
3)确定其他类型约束条件的系数矩阵及常数向量
 如本例中的优化模型所示,容易确定其余的输入参数,线性不等式约束条件的系数矩阵A和常数向量分别为: A=[1 1],b=[2 ],线性等式约束不存在,故Aeq=[],beq=[],设计变量X的上、下界向量:lb=[0 0]',ub=[inf inf]',其中inf表示无穷大。
4)调用fmincon函数进行求解
经过上述各步骤设置以后,可以编制主程序进行优化求解,相应的代码如下:
>> x0=[1 1]; %设置计算初始值
>> options=optimset('LargeScale','off','display','iter'); %设定优化选项参数
>> [x,fval,exitflag]=fmincon(@myobj,x0,A,b,[],[],lb,ub,@mycon,options) %进行优化求解
讲过运算以后得到结果如下所示:
Optimization terminated successfully:
 First-order optimality measure less than options.TolFun and
  maximum constraint violation is less than options.TolCon
Active Constraints:
     3
     4
x =
    1.1190    0.8810
fval =
   -7.6771
exitflag =
     1
Read more ...

MATLAB commands and functions list

A a

abs 绝对值、模、字符的ASCII码值
acos 反余弦
acosh 反双曲余弦
acot 反余切
acoth 反双曲余切
acsc 反余割
acsch 反双曲余割
align 启动图形对象几何位置排列工具
all 所有元素非零为真
angle 相角
ans 表达式计算结果的缺省变量名
any 所有元素非全零为真
area 面域图
argnames 函数M文件宗量名
asec 反正割
asech 反双曲正割
asin 反正弦
asinh 反双曲正弦
assignin 向变量赋值
atan 反正切
atan2 四象限反正切
atanh 反双曲正切
autumn 红黄调秋色图阵
axes 创建轴对象的低层指令
axis 控制轴刻度和风格的高层指令


B b

bar 二维直方图
bar3 三维直方图
bar3h 三维水平直方图
barh 二维水平直方图
base2dec X进制转换为十进制
bin2dec 二进制转换为十进制
blanks 创建空格串
bone 蓝色调黑白色图阵
box 框状坐标轴
break while 或for 环中断指令
brighten 亮度控制


C c

capture (3版以前)捕获当前图形
cart2pol 直角坐标变为极或柱坐标
cart2sph 直角坐标变为球坐标
cat 串接成高维数组
caxis 色标尺刻度
cd 指定当前目录
cdedit 启动用户菜单、控件回调函数设计工具
cdf2rdf 复数特征值对角阵转为实数块对角阵
ceil 向正无穷取整
cell 创建元胞数组
cell2struct 元胞数组转换为构架数组
celldisp 显示元胞数组内容
cellplot 元胞数组内部结构图示
char 把数值、符号、内联类转换为字符对象
chi2cdf 分布累计概率函数
chi2inv 分布逆累计概率函数
chi2pdf 分布概率密度函数
chi2rnd 分布随机数发生器
chol Cholesky分解
clabel 等位线标识
cla 清除当前轴
class 获知对象类别或创建对象
clc 清除指令窗
clear 清除内存变量和函数
clf 清除图对象
clock 时钟
colorcube 三浓淡多彩交叉色图矩阵
colordef 设置色彩缺省值
colormap 色图
colspace 列空间的基
close 关闭指定窗口
colperm 列排序置换向量
comet 彗星状轨迹图
comet3 三维彗星轨迹图
compass 射线图
compose 求复合函数
cond (逆)条件数
condeig 计算特征值、特征向量同时给出条件数
condest 范 -1条件数估计
conj 复数共轭
contour 等位线
contourf 填色等位线
contour3 三维等位线
contourslice 四维切片等位线图
conv 多项式乘、卷积
cool 青紫调冷色图
copper 古铜调色图
cos 余弦
cosh 双曲余弦
cot 余切
coth 双曲余切
cplxpair 复数共轭成对排列
csc 余割
csch 双曲余割
cumsum 元素累计和
cumtrapz 累计梯形积分
cylinder 创建圆柱


D d

dblquad 二重数值积分
deal 分配宗量
deblank 删去串尾部的空格符
dec2base 十进制转换为X进制
dec2bin 十进制转换为二进制
dec2hex 十进制转换为十六进制
deconv 多项式除、解卷
delaunay Delaunay 三角剖分
del2 离散Laplacian差分
demo Matlab演示
det 行列式
diag 矩阵对角元素提取、创建对角阵
diary Matlab指令窗文本内容记录
diff 数值差分、符号微分
digits 符号计算中设置符号数值的精度
dir 目录列表
disp 显示数组
display 显示对象内容的重载函数
dlinmod 离散系统的线性化模型
dmperm 矩阵Dulmage-Mendelsohn 分解
dos 执行DOS 指令并返回结果
double 把其他类型对象转换为双精度数值
drawnow 更新事件队列强迫Matlab刷新屏幕
dsolve 符号计算解微分方程


E e

echo M文件被执行指令的显示
edit 启动M文件编辑器
eig 求特征值和特征向量
eigs 求指定的几个特征值
end 控制流FOR等结构体的结尾元素下标
eps 浮点相对精度
error 显示出错信息并中断执行
errortrap 错误发生后程序是否继续执行的控制
erf 误差函数
erfc 误差补函数
erfcx 刻度误差补函数
erfinv 逆误差函数
errorbar 带误差限的曲线图
etreeplot 画消去树
eval 串演算指令
evalin 跨空间串演算指令
exist 检查变量或函数是否已定义
exit 退出Matlab环境
exp 指数函数
expand 符号计算中的展开操作
expint 指数积分函数
expm 常用矩阵指数函数
expm1 Pade法求矩阵指数
expm2 Taylor法求矩阵指数
expm3 特征值分解法求矩阵指数
eye 单位阵
ezcontour 画等位线的简捷指令
ezcontourf 画填色等位线的简捷指令
ezgraph3 画表面图的通用简捷指令
ezmesh 画网线图的简捷指令
ezmeshc 画带等位线的网线图的简捷指令
ezplot 画二维曲线的简捷指令
ezplot3 画三维曲线的简捷指令
ezpolar 画极坐标图的简捷指令
ezsurf 画表面图的简捷指令
ezsurfc 画带等位线的表面图的简捷指令



F f

factor 符号计算的因式分解
feather 羽毛图
feedback 反馈连接
feval 执行由串指定的函数
fft 离散Fourier变换
fft2 二维离散Fourier变换
fftn 高维离散Fourier变换
fftshift 直流分量对中的谱
fieldnames 构架域名
figure 创建图形窗
fill3 三维多边形填色图
find 寻找非零元素下标
findobj 寻找具有指定属性的对象图柄
findstr 寻找短串的起始字符下标
findsym 机器确定内存中的符号变量
finverse 符号计算中求反函数
fix 向零取整
flag 红白蓝黑交错色图阵
fliplr 矩阵的左右翻转
flipud 矩阵的上下翻转
flipdim 矩阵沿指定维翻转
floor 向负无穷取整
flops 浮点运算次数
flow Matlab提供的演示数据
fmin 求单变量非线性函数极小值点(旧版)
fminbnd 求单变量非线性函数极小值点
fmins 单纯形法求多变量函数极小值点(旧版)
fminunc 拟牛顿法求多变量函数极小值点
fminsearch 单纯形法求多变量函数极小值点
fnder 对样条函数求导
fnint 利用样条函数求积分
fnval 计算样条函数区间内任意一点的值
fnplt 绘制样条函数图形
fopen 打开外部文件
for 构成for环用
format 设置输出格式
fourier Fourier 变换
fplot 返函绘图指令
fprintf 设置显示格式
fread 从文件读二进制数据
fsolve 求多元函数的零点
full 把稀疏矩阵转换为非稀疏阵
funm 计算一般矩阵函数
funtool 函数计算器图形用户界面
fzero 求单变量非线性函数的零点


G g

gamma 函数
gammainc 不完全 函数
gammaln 函数的对数
gca 获得当前轴句柄
gcbo 获得正执行"回调"的对象句柄
gcf 获得当前图对象句柄
gco 获得当前对象句柄
geomean 几何平均值
get 获知对象属性
getfield 获知构架数组的域
getframe 获取影片的帧画面
ginput 从图形窗获取数据
global 定义全局变量
gplot 依图论法则画图
gradient 近似梯度
gray 黑白灰度
grid 画分格线
griddata 规则化数据和曲面拟合
gtext 由鼠标放置注释文字
guide 启动图形用户界面交互设计工具


H h

harmmean 调和平均值
help 在线帮助
helpwin 交互式在线帮助
helpdesk 打开超文本形式用户指南
hex2dec 十六进制转换为十进制
hex2num 十六进制转换为浮点数
hidden 透视和消隐开关
hilb Hilbert矩阵
hist 频数计算或频数直方图
histc 端点定位频数直方图
histfit 带正态拟合的频数直方图
hold 当前图上重画的切换开关
horner 分解成嵌套形式
hot 黑红黄白色图
hsv 饱和色图


I i

if-else-elseif 条件分支结构
ifft 离散Fourier反变换
ifft2 二维离散Fourier反变换
ifftn 高维离散Fourier反变换
ifftshift 直流分量对中的谱的反操作
ifourier Fourier反变换
i, j 缺省的"虚单元"变量
ilaplace Laplace反变换
imag 复数虚部
image 显示图象
imagesc 显示亮度图象
imfinfo 获取图形文件信息
imread 从文件读取图象
imwrite 把图象写成文件
ind2sub 单下标转变为多下标
inf 无穷大
info MathWorks公司网点地址
inline 构造内联函数对象
inmem 列出内存中的函数名
input 提示用户输入
inputname 输入宗量名
int 符号积分
int2str 把整数数组转换为串数组
interp1 一维插值
interp2 二维插值
interp3 三维插值
interpn N维插值
interpft 利用FFT插值
intro Matlab自带的入门引导
inv 求矩阵逆
invhilb Hilbert矩阵的准确逆
ipermute 广义反转置
isa 检测是否给定类的对象
ischar 若是字符串则为真
isequal 若两数组相同则为真
isempty 若是空阵则为真
isfinite 若全部元素都有限则为真
isfield 若是构架域则为真
isglobal 若是全局变量则为真
ishandle 若是图形句柄则为真
ishold 若当前图形处于保留状态则为真
isieee 若计算机执行IEEE规则则为真
isinf 若是无穷数据则为真
isletter 若是英文字母则为真
islogical 若是逻辑数组则为真
ismember 检查是否属于指定集
isnan 若是非数则为真
isnumeric 若是数值数组则为真
isobject 若是对象则为真
isprime 若是质数则为真
isreal 若是实数则为真
isspace 若是空格则为真
issparse 若是稀疏矩阵则为真
isstruct 若是构架则为真
isstudent 若是Matlab学生版则为真
iztrans 符号计算Z反变换


J j , K k

jacobian 符号计算中求Jacobian 矩阵
jet 蓝头红尾饱和色
jordan 符号计算中获得 Jordan标准型
keyboard 键盘获得控制权
kron Kronecker乘法规则产生的数组


L l

laplace Laplace变换
lasterr 显示最新出错信息
lastwarn 显示最新警告信息
leastsq 解非线性最小二乘问题(旧版)
legend 图形图例
lighting 照明模式
line 创建线对象
lines 采用plot 画线色
linmod 获连续系统的线性化模型
linmod2 获连续系统的线性化精良模型
linspace 线性等分向量
ln 矩阵自然对数
load 从MAT文件读取变量
log 自然对数
log10 常用对数
log2 底为2的对数
loglog 双对数刻度图形
logm 矩阵对数
logspace 对数分度向量
lookfor 按关键字搜索M文件
lower 转换为小写字母
lsqnonlin 解非线性最小二乘问题
lu LU分解


M m

mad 平均绝对值偏差
magic 魔方阵
maple &nb, sp; 运作 Maple格式指令
mat2str 把数值数组转换成输入形态串数组
material 材料反射模式
max 找向量中最大元素
mbuild 产生EXE文件编译环境的预设置指令
mcc 创建MEX或EXE文件的编译指令
mean 求向量元素的平均值
median 求中位数
menuedit 启动设计用户菜单的交互式编辑工具
mesh 网线图
meshz 垂帘网线图
meshgrid 产生"格点"矩阵
methods 获知对指定类定义的所有方法函数
mex 产生MEX文件编译环境的预设置指令
mfunlis 能被mfun计算的MAPLE经典函数列表
mhelp 引出 Maple的在线帮助
min 找向量中最小元素
mkdir 创建目录
mkpp 逐段多项式数据的明晰化
mod 模运算
more 指令窗中内容的分页显示
movie 放映影片动画
moviein 影片帧画面的内存预置
mtaylor 符号计算多变量Taylor级数展开


N n

ndims 求数组维数
NaN 非数(预定义)变量
nargchk 输入宗量数验证
nargin 函数输入宗量数
nargout 函数输出宗量数
ndgrid 产生高维格点矩阵
newplot 准备新的缺省图、轴
nextpow2 取最接近的较大2次幂
nnz 矩阵的非零元素总数
nonzeros 矩阵的非零元素
norm 矩阵或向量范数
normcdf 正态分布累计概率密度函数
normest 估计矩阵2范数
norminv 正态分布逆累计概率密度函数
normpdf 正态分布概率密度函数
normrnd 正态随机数发生器
notebook 启动Matlab和Word的集成环境
null 零空间
num2str 把非整数数组转换为串
numden 获取最小公分母和相应的分子表达式
nzmax 指定存放非零元素所需内存


O o

ode1 非Stiff 微分方程变步长解算器
ode15s Stiff 微分方程变步长解算器
ode23t 适度Stiff 微分方程解算器
ode23tb Stiff 微分方程解算器
ode45 非Stiff 微分方程变步长解算器
odefile ODE 文件模板
odeget 获知ODE 选项设置参数
odephas2 ODE 输出函数的二维相平面图
odephas3 ODE 输出函数的三维相空间图
odeplot ODE 输出函数的时间轨迹图
odeprint 在Matlab指令窗显示结果
odeset 创建或改写 ODE选项构架参数值
ones 全1数组
optimset 创建或改写优化泛函指令的选项参数值
orient 设定图形的排放方式
orth 值空间正交化


P p

pack 收集Matlab内存碎块扩大内存
pagedlg 调出图形排版对话框
patch 创建块对象
path 设置Matlab搜索路径的指令
pathtool 搜索路径管理器
pause 暂停
pcode 创建预解译P码文件
pcolor 伪彩图
peaks Matlab提供的典型三维曲面
permute 广义转置
pi (预定义变量)圆周率
pie 二维饼图
pie3 三维饼图
pink 粉红色图矩阵
pinv 伪逆
plot 平面线图
plot3 三维线图
plotmatrix 矩阵的散点图
plotyy 双纵坐标图
poissinv 泊松分布逆累计概率分布函数
poissrnd 泊松分布随机数发生器
pol2cart 极或柱坐标变为直角坐标
polar 极坐标图
poly 矩阵的特征多项式、根集对应的多项式
poly2str 以习惯方式显示多项式
poly2sym 双精度多项式系数转变为向量符号多项式
polyder 多项式导数
polyfit 数据的多项式拟合
polyval 计算多项式的值
polyvalm 计算矩阵多项式
pow2 2的幂
ppval 计算分段多项式
pretty 以习惯方式显示符号表达式
print 打印图形或SIMULINK模型
printsys 以习惯方式显示有理分式
prism 光谱色图矩阵
procread 向MAPLE输送计算程序
profile 函数文件性能评估器
propedit 图形对象属性编辑器
pwd 显示当前工作目录


Q q

quad 低阶法计算数值积分
quad8 高阶法计算数值积分(QUADL)
quit 推出Matlab 环境
quiver 二维方向箭头图
quiver3 三维方向箭头图


R r

rand 产生均匀分布随机数
randn 产生正态分布随机数
randperm 随机置换向量
range 样本极差
rank 矩阵的秩
rats 有理输出
rcond 矩阵倒条件数估计
real 复数的实部
reallog 在实数域内计算自然对数
realpow 在实数域内计算乘方
realsqrt 在实数域内计算平方根
realmax 最大正浮点数
realmin 最小正浮点数
rectangle 画"长方框"
rem 求余数
repmat 铺放模块数组
reshape 改变数组维数、大小
residue 部分分式展开
return 返回
ribbon 把二维曲线画成三维彩带图
rmfield 删去构架的域
roots 求多项式的根
rose 数扇形图
rot90 矩阵旋转90度
rotate 指定的原点和方向旋转
rotate3d 启动三维图形视角的交互设置功能
round 向最近整数圆整
rref 简化矩阵为梯形形式
rsf2csf 实数块对角阵转为复数特征值对角阵
rsums Riemann和


S s

save 把内存变量保存为文件
scatter 散点图
scatter3 三维散点图
sec 正割
sech 双曲正割
semilogx X轴对数刻度坐标图
semilogy Y轴对数刻度坐标图
series 串联连接
set 设置图形对象属性
setfield 设置构架数组的域
setstr 将ASCII码转换为字符的旧版指令
sign 根据符号取值函数
signum 符号计算中的符号取值函数
sim 运行SIMULINK模型
simget 获取SIMULINK模型设置的仿真参数
simple 寻找最短形式的符号解
simplify 符号计算中进行简化操作
simset 对SIMULINK模型的仿真参数进行设置
simulink 启动SIMULINK模块库浏览器
sin 正弦
sinh 双曲正弦
size 矩阵的大小
slice 立体切片图
solve 求代数方程的符号解
spalloc 为非零元素配置内存
sparse 创建稀疏矩阵
spconvert 把外部数据转换为稀疏矩阵
spdiags 稀疏对角阵
spfun 求非零元素的函数值
sph2cart 球坐标变为直角坐标
sphere 产生球面
spinmap 色图彩色的周期变化
spline 样条插值
spones 用1置换非零元素
sprandsym 稀疏随机对称阵
sprank 结构秩
spring 紫黄调春色图
sprintf 把格式数据写成串
spy 画稀疏结构图
sqrt 平方根
sqrtm 方根矩阵
squeeze 删去大小为1的"孤维"
sscanf 按指定格式读串
stairs 阶梯图
std 标准差
stem 二维杆图
step 阶跃响应指令
str2double 串转换为双精度值
str2mat 创建多行串数组
str2num 串转换为数
strcat 接成长串
strcmp 串比较
strjust 串对齐
strmatch 搜索指定串
strncmp 串中前若干字符比较
strrep 串替换
strtok 寻找第一间隔符前的内容
struct 创建构架数组
struct2cell 把构架转换为元胞数组
strvcat 创建多行串数组
sub2ind 多下标转换为单下标
subexpr 通过子表达式重写符号对象
subplot 创建子图
subs 符号计算中的符号变量置换
subspace 两子空间夹角
sum 元素和
summer 绿黄调夏色图
superiorto 设定优先级
surf 三维着色表面图
surface 创建面对象
surfc 带等位线的表面图
surfl 带光照的三维表面图
surfnorm 空间表面的法线
svd 奇异值分解
svds 求指定的若干奇异值
switch-case-otherwise 多分支结构
sym2poly 符号多项式转变为双精度多项式系数向量
symmmd 对称最小度排序
symrcm 反向Cuthill-McKee排序
syms 创建多个符号对象


T t

tan 正切
tanh 双曲正切
taylortool 进行Taylor逼近分析的交互界面
text 文字注释
tf 创建传递函数对象
tic 启动计时器
title 图名
toc 关闭计时器
trapz 梯形法数值积分
treelayout 展开树、林
treeplot 画树图
tril 下三角阵
trim 求系统平衡点
trimesh 不规则格点网线图
trisurf 不规则格点表面图 triu 上三角阵 try-catch 控制流中的Try-catch结构 type 显示M文件
U u
uicontextmenu 创建现场菜单
uicontrol 创建用户控件
uimenu 创建用户菜单
unmkpp 逐段多项式数据的反明晰化
unwrap 自然态相角
upper 转换为大写字母


V v

var 方差
varargin 变长度输入宗量
varargout 变长度输出宗量
vectorize 使串表达式或内联函数适于数组运算
ver 版本信息的获取
view 三维图形的视角控制
voronoi Voronoi多边形
vpa 任意精度(符号类)数值


W w

warning 显示警告信息
what 列出当前目录上的文件
whatsnew 显示Matlab中 Readme文件的内容
which 确定函数、文件的位置
while 控制流中的While环结构
white 全白色图矩阵
whitebg 指定轴的背景色
who 列出内存中的变量名
whos 列出内存中变量的详细信息
winter 蓝绿调冬色图
workspace 启动内存浏览器


X x , Y y , Z z

xlabel X轴名
xor 或非逻辑
yesinput 智能输入指令
ylabel Y轴名
zeros 全零数组
zlabel Z轴名
zoom 图形的变焦放大和缩小
ztrans 符号计算Z变换
Read more ...

MATLAB 矩阵相关 command

pascal(m,n)帕斯卡矩阵
fliplr(a) 矩阵左右翻转
flipud(a) 矩阵上下翻转
rot90(a)
rot90(a,k) 矩阵逆时针旋转90度(把你的头顺时针旋转90看原数就可以知道结果了,^-^)
k参数定义为逆时针旋转90*k度。
flipdim(a,k) 矩阵对应维数数值翻转,如k=1时,行(上下)翻转,k=2时,列(左右)翻转。
tril(a)
tril(a,k) 矩阵的下三角部分(包括对角线元素),对应k=0时的取值数。
k参数设置为正负数值对应对角线向上或向下移动k行划分下三角元素。
triu(a)
tril(a,k) 矩阵的上三角部分(包括对角线元素),对应k=0时的取值数。
k参数设置为正负数值对应对角线向上或向下移动k行划分上三角元素。
diag(a)
diag(a,k) 生成对角矩阵或取出对角元素,对应k=0时的取值数。
k参数设置为正负数值对应对角线向上或向下移动k行取对角元素或生成对角矩阵。
repmat(a,m,n) 矩阵复制,把矩阵a作为一个单位计算,复制成m*n的矩阵,其每一元素都含一个矩阵a,实际结果为一个size(a,1)*m行,size(a,2)*n列的矩阵。
w=meshgrid(s,t)
[u,v]=meshgrid(s,t) 生成行m=size(t,1)*size(t,2),列n=size(s,1)*size(s,2))阶的两个矩阵。其中u为按行顺序取s的n个矩阵元数,按列排列重复m行,v为按列顺序取t的m个矩阵元数 ,按行排列重复n列。只生成一个矩阵时,w=u。
eye(a)
eye(a,k) 生成a阶单位方阵
k参数设置为生成a×k阶单位矩阵,即生成a阶单位方阵后,取前k列,不足补0。
ones(a)
ones(a,k) 生成a阶全1方阵
k参数设置生成a×k阶全1矩阵。
zeros(a)
zeros(a,k) 生成a阶全0方阵
k参数设置生成a×k阶全0矩阵。
inv(a) 生成a的逆矩阵
Read more ...

Linux 下安装 MATLAB

$ sudo mount -o loop /home/username/Downloads/matlabR2012b.iso /mnt
$ cd /mnt
$ sudo ./install

run Matlab
$ matlab -desktop
or
$ /usr/local/MATLAB/R2012b/bin/matlab -desktop

解决 /lib/libc.so.6: not found

For 64 bit:
$ sudo ln -s /lib/x86_64-linux-gnu/libc-2.13.so /lib64/libc.so.6 

For 32 bit:
$ sudo ln -s /lib/i386-linux-gnu/libc-2.13.so /lib/libc.so.6
Read more ...

How do I extract data from MATLAB figures?

Subject:

How do I extract data from MATLAB figures?

Problem Description:

I have a few MATLAB figures, but no MATLAB code associated with it. I want to extract the data from the curves in the figures.

Solution:

Below is a step by step example to extract data from the curve in a MATLAB figure :

Assume that the figure is stored in a file called 'example.fig'.

<step 1>
Open the figure file:

open('example.fig');

%or
figure;
plot(1:10)


<step 2>
Get a handle to the current figure:

h = gcf;


<step 3>
The data that is plotted is usually a 'child' of the Axes object. The axes objects are themselves children of the figure. You can go down their hierarchy as follows:

axesObjs = get(h, 'Children');
dataObjs = get(axesObjs, 'Children'); 


<step 4>
Extract values from the dataObjs of your choice. You can check their type by typing:

objTypes = get(dataObjs, 'Type');


* NOTE : Different objects like 'Line' and 'Surface' will store data differently. Based on the 'Type', you can search the documentation for how each type stores its data.

<step 5>
Lines of code similar to the following would be required to bring the data to MATLAB Workspace:
xdata = get(dataObjs, 'XData');
ydata = get(dataObjs, 'YData');
zdata = get(dataObjs, 'ZData');
Read more ...

数理统计与Matlab: 第5章 方差分析

This post is from here.


5.1 单因素方差分析

5.1.1 方差分析的基本概念
在实际问题中,人们常常需要在不同的条件下对所研究的对象进行对比试验,从而得到若干组数据(样本)。方差分析就是一种分析、处理多组实验数据间均值差异的显著性的统计方法。其主要任务是,通过对数据的分析处理,搞清楚各实验条件对实验结果的影响,以便更有效地指导实践,提高经济效益或者科研水平。
在统计中,人们称受控制的条件为因素,因素所处的状态称为水平
如果只让一个因素变动,取该因素的多个不同水平进行试验,而其他因素保持不变,称该试验为单因素试验。例如小麦种植产量,只考虑"品种"这一因素,研究4个不同品种产量的差异,其它诸如施肥方案、灌溉方案等因素保持一致,就是一个4水平单因素试验。
如果同时考虑两个因素,例如4个小麦品种在3种不同施肥方案下的产量,就是一个双因素试验。
对于组实验数据,我们假定都来自正态总体,并且具有相同的方差(称为方差齐性),要检验这相互独立的个正态总体
 
均值间有无差异,即:
H0; H1:诸不全相同
前面我们讲过两正态总体均值的假设检验,有T检验的方法。自然有一个想法,对于,分别检验是否成立,若所有搭配均不拒绝,则接受H0,只要有一种搭配拒绝原假设认为,那就拒绝H0,看起来也不算麻烦。不妨称上述想法为"两两T检验法"。
回忆前面内容,设为来自正态总体的简单随机样本,为来自正态总体的简单随机样本,且两样本独立。为比较两个总体的期望,提出如下原假设:
H0
H0成立时,检验统计量
我们给出函数t2test.m,解决上述计算问题。
function T=t2test(x,y)
m=length(x);
n=length(y);
vx=var(x);
vy=var(y);
a=(mean(x)-mean(y));
b=m+n-2;
c=(m-1)*vx+(n-1)*vy;
d=sqrt(m*n/(m+n));
T=a*d*sqrt(b/c);
以下给出m=10,n=10,且两总体皆服从标准正态分布的情形下,万次模拟的拒绝频率。以下命令文件保存为PnT2.m
N=10000;
m=10; n=10;
alpha=0.05;
t0=tinv(1-alpha/2,m+n-2);
P=0;
for k=1:N
x=randn(1,10); y=randn(1,10);
T=t2test(x,y);
if abs(T)>t0
P=P+1;
end
end
P=P/N
执行上述程序,发现每次频率都在0.05附近,说明上述两个正态总体均值的T检验的确是水平为的检验。
我们设想有8组数据,客观上都是来自标准正态分布,没有差异,每组样本容量都是10。现在用前述"两两T检验法"进行检验,下述程序计算出了万次模拟中拒绝的频率。
N=10000;
n=10; r=8;
alpha=0.05;
t0=tinv(1-alpha/2,n+n-2);
P=0;
for k=1:N
x=randn(8,10);
E=mean(x,2);
[EE,I]=sort(E);
X=x(I,:);
T=t2test(X(1,:),X(8,:));
if abs(T)>t0
P=P+1;
end
end
P=P/N
上述程序模拟发现,拒绝频率大约在0.45左右,严重偏离0.05,说明依照"两两T检验"犯第一类错误的概率严重增大,判定结果很不可靠。
对于8组数据,两两比较共种组合,若每种组合接受原假设的概率为0.95,则28种组合都接受原假设的概率大致估计为,拒绝概率大致估计为0.76。由于相关性,拒绝概率没有达到0.76,但0.45也相当大了。
为了避免上述问题的出现,1923年,波兰数学家R.A.Fisher提出了方差分析(Analysis of Variance简称ANOVA) 法,可以同时判定多组数据均值间差异的显著性检验问题。其检验统计量在H0成立时服从F分布,这里F分布就是以Fisher姓氏的第一个字母命名的。
5.1.2 单因素方差分析的计算
设有组数据,表示因素A的个水平,每组有个观测值。我们已知实际结果具有以下结构:
 ()
表示水平Ai下的理论均值,为实验误差,诸相互独立且服从正态分布
为了看出因素A个水平影响的大小,将进行分解,令
, 
表示水平Ai对试验结果的影响,称为Ai的水平效应。显然
这时数据有如下结构:
 () (5-1)
于是,我们需要进行的假设检验为:
H0; H1:诸不全为零 (5-2)

,() , 
 (5-3)
总离差平方和,它反映了样本观测值之间的总的变异程度。以下我们将分解为两部分,以便区别水平效应与随机误差的影响。
其中
,  (5-4)
组内平方和,它反映了每组的组内随机误差。称组间平方和,反应的是组与组之间的差异。上述推导说明,总离差平方可以分解为
 (5-5)
一个自然的想法是:如果在总离差平方和中,所占比例很大,则拒绝原假设,认为客观上存在水平效应。
H0成立时容易计算
 ;  (5-6)
因此,当H0成立时,有
,  (5-7)
 (5-8)
对于自由度,求临界值,当时拒绝H0即可。
表5-1 单因素方差分析表
方差来源
平方和
自由度
均方
F值
临界值
显著性
组间
SA
 
误差
SE
总和
ST
 
实际计算时常采用方差分析表,如表5-1所示。当时,称为不显著,即认为各组均值之间没有显著差异,在显著性一栏不做任何标记。当时,称为较显著,即认为各组均值之间有较显著差异,在显著性一栏用(*)标记。当时,称为显著,即认为各组均值之间有显著差异,在显著性一栏用*标记。当时,称为极显著,即认为各组均值之间有极显著差异,在显著性一栏用**标记。
上述传统的方差计算表,在计算机普及后稍有变动,表中最后两列可以变动为直接计算H0成立时F分布大于此F值的概率,是否显著一看自明。
例5.1 为了研究三种不同伤寒杆菌对于小白鼠存活天数的影响,分三组实验,实验数据如下:
A1: 2, 4, 3, 2, 4, 7, 7, 2, 5, 4
A2: 5, 6, 8, 5, 10,7, 12,6, 6
A3: 7, 11,6, 6, 7, 9, 5, 10,6, 3,10
试检验不同伤寒杆菌对于小白鼠存活天数有无显著影响?
 原假设H0没有显著差异。以下利用Matlab计算,并将计算结果汇总到表5-2中。
x1=[2, 4, 3, 2, 4, 7, 7, 2, 5, 4];
x2=[5, 6, 8, 5, 10,7, 12,6, 6];
x3=[7, 11,6, 6, 7, 9, 5, 10,6, 3,10];
n1=length(x1);
n2=length(x2);
n3=length(x3);
X1=sum(x1), mx1=mean(x1),
X2=sum(x2), mx2=mean(x2),
X3=sum(x3), mx3=mean(x3),
n=n1+n2+n3, X=X1+X2+X3, mx=X/n,
SE=(x1-mx1)*(x1-mx1)'+(x2-mx2)*(x2-mx2)'+(x3-mx3)*(x3-mx3)',
SA=n1*(mx1-mx)^2+n2*(mx2-mx)^2+n3*(mx3-mx)^2,
F=(SA/2)/(SE/27),
finv(0.9,2,27),
finv(0.95,2,27),
finv(0.99,2,27),
p=1-fcdf(F,2,27),
表5-2 不同伤寒杆菌对小白鼠存活天数影响的方差分析表
方差来源
平方和
自由度
均方
F值
临界值
显著性
组间
70.4293
2
35.2146
6.9030
2.5106
**
误差
137.7374
27
5.1014
3.3541
总和
208.1667
29
 
5.4881
可以认为不同伤寒杆菌对小白鼠存活天数有极显著影响。
上述最后一行命令执行的结果为p=0.0038,实际判定时,此值更能说明显著性程度。
当各组样本容量相等时,Matlab自带了单因素方差分析函数anova1,并且自动返回类似的方差分析表,调用格式为anova1(x),这里x为矩阵,按列分组,第一列为第一组数据,要求每组数据相同。
例5.2 某班学生共分别住在四个宿舍,某次英语水平考试成绩如下:
表5-3 各宿舍学生英语成绩表
宿舍一
86
74
78
76
73
82
宿舍二
59
76
63
66
75
65
宿舍三
79
71
82
66
87
96
宿舍四
86
91
95
77
88
85
问不同宿舍学生英语水平有无显著差异?
 利用复制粘贴的办法输入矩阵x,由于anova1要求按列分组,故使用x=x'命令,然后再执行命令
p=anova1(x),
返回值为p =0.002,故差异极显著。Matlab同时返回了两个图形。
图5-1 Matlab中anova1返回的方差分析表

图5-2 Matlab中anova1返回的各组数据图
图5-2中各线含义依次为:最小值、1/4分位数、中位数、3/4分位数、最大值。
5.1.3 单因素方差分析的多重比较
经过方差分析之后,如果拒绝原假设,认为各组之间的均值有显著差异,那么,这个判断是对整体而言的,并不是说每两个不同的组之间均值都存在显著差异。那么,如何确定哪两个组之间有显著差异、无显著差异呢?这就要对每种搭配做一对一的比较,即多重比较。Matlab提供了multcompare函数用于进行多重比较,为了使用这个函数,在用anova1做方差分析的时候要使用如下的三个输出的调用格式
[P,ANOVATAB,STATS] = anova1(x)
在例5.2的数据输入之后,假定x已经经过转置使得每组数据按列排列,执行上述命令后,则除了返回上述两个图形外,还返回
P =
0.0020
ANOVATAB =
'Source' 'SS' 'df' 'MS' 'F' 'Prob>F'
'Columns' [1.1963e+003] [ 3] [398.7778] [7.0768] [0.0020]
'Error' [1.1270e+003] [20] [ 56.3500] [] []
'Total' [2.3233e+003] [23] [] [] []
STATS =
gnames: [4x1 char]
n: [6 6 6 6]
source: 'anova1'
means: [78.1667 67.3333 80.1667 87]
df: 20
s: 7.5067
其中最后一个输出结果STATS用于多重比较函数调用。继续执行命令
COMPARISON = multcompare(STATS,'alpha',0.1)
得到输出结果

COMPARISON =
1.0000 2.0000 0.2251 10.8333 21.4415
1.0000 3.0000 -12.6082 -2.0000 8.6082
1.0000 4.0000 -19.4415 -8.8333 1.7749
2.0000 3.0000 -23.4415 -12.8333 -2.2251
2.0000 4.0000 -30.2749 -19.6667 -9.0585
3.0000 4.0000 -17.4415 -6.8333 3.7749
上述结果中,逐行显示一对一比较的结果。第一列与第二列是参与比较的组号,第三列与第五列是均值差的置信区间,置信水平由输入变量alpha=0.1确定。第四列为两组样本均值的差。若第三列置信下限与第五列置信上限正负符号相同,则差异显著。由于alpha=0.1,我们可以发现第1,2组差异较显著,第2,3组差异较显著,第2,4组差异较显著,其它对比差异不显著。
再修改上述输入变量alpha的值,重新计算
COMPARISON = multcompare(STATS,'alpha',0.05)
输出结果为
COMPARISON =
1.0000 2.0000 -1.2972 10.8333 22.9639
1.0000 3.0000 -14.1305 -2.0000 10.1305
1.0000 4.0000 -20.9639 -8.8333 3.2972
2.0000 3.0000 -24.9639 -12.8333 -0.7028
2.0000 4.0000 -31.7972 -19.6667 -7.5361
3.0000 4.0000 -18.9639 -6.8333 5.2972
只有第2,3组,第2,4组有显著差异。继续执行
COMPARISON = multcompare(STATS,'alpha',0.01)
输出结果为
COMPARISON =
1.0000 2.0000 -4.5448 10.8333 26.2115
1.0000 3.0000 -17.3781 -2.0000 13.3781
1.0000 4.0000 -24.2115 -8.8333 6.5448
2.0000 3.0000 -28.2115 -12.8333 2.5448
2.0000 4.0000 -35.0448 -19.6667 -4.2885
3.0000 4.0000 -22.2115 -6.8333 8.5448
发现只有第2,4组有极显著的差异。

5.2 双因素方差分析

5.2.1 有重复实验的双因素方差分析
如果实验同时受到两个因素的共同影响,因素个水平,因素个水平,一共有种搭配方案,每种搭配有做个实验数据,实验数据表依照表5-4所示排列。
首先引进记号
的三个下标中,若有的下标为点,则表示对该下标求和,例如
, , 
相应的平均值记为
, , 
依此类推。
将上述诸看成是随机变量,并且假设满足以下线性模型
,  (5-9)
表5-4 双因素有重复试验数据表
 
因素
 
 
其中是各种搭配下的总平均数的理论值,是因素的第个水平的效应,是因素的第个水平的效应,的交互作用,即搭配效应。是相互独立的正态随机变量,均值为零,方差为,并且诸相互独立。满足
 (5-10)
对于上述有重复试验的双因素实验,方差分析检验如下三个假设:
H01 (5-11)
H02 (5-12)
H03 (5-13)
类似单因素方差分析,记总离差平方和为
 (5-14)
则有如下分解式
 自由度为 (5-15)
其中
, 自由度为 (5-16)
, 自由度为 (5-17)
, 自由度为 (5-18)
, 自由度为 (5-19)
由此可得检验统计量,当H01成立时,
 (5-20)
时,拒绝H01,认为因素作用显著。
H02成立时,
 (5-21)
时,拒绝H02,认为因素作用显著。
H03成立时,
 (5-22)
时,拒绝H03,认为交互作用显著。
方差计算可以类似地用方差计算表进行。
Matlab自带的函数anova2用于处理双因素方差分析,调用格式为
[P,TABLE,STATS] = anova2(x,n)
其中P返回概率值,TABLE返回方差计算表,STATS返回的信息用于多重比较。输入变量n表示每种搭配下样本容量,记号同前述公式。x为数据矩阵,为na行b列,格式与表5-4完全一致。
例5.3 杨树一年中生长高度受两种因素影响,A:施肥方案,B:深翻方案。对于4种施肥方案及3种深翻方案,共12种搭配,各实验3株,实验结果如表5-5所示。试问施肥方案、深翻方案、两者的交互作用对于苗高有无显著影响?
表5-5 杨树苗增高实验数据表
 
B1
B2
B3
A1
52 43 39
41 47 53
49 38 42
A2
48 37 29
50 41 30
36 48 47
A3
34 42 38
36 39 44
37 40 32
A4
45 58 42
44 46 60
43 56 41

 利用复制粘贴的办法对矩阵x赋值,注意到anova2的要求,重复数据按列排列,故对x进行转置,x=x'。转置之后x为9行4列矩阵,每列表示4种施肥方案的数据,行对应的是3种深翻方案。执行Matlab命令:
[P,TABLE,STATS] = anova2(x,3)
计算结果为:
P =
0.0259 0.7504 0.9540
TABLE =
'Source' 'SS' 'df' 'MS' 'F' 'Prob>F'
'Columns' [ 562.0833] [ 3] [187.3611] [3.6838] [0.0259]
'Rows' [ 29.5556] [ 2] [ 14.7778] [0.2906] [0.7504]
'Interaction' [ 76.6667] [ 6] [ 12.7778] [0.2512] [0.9540]
'Error' [1.2207e+003] [24] [ 50.8611] [] []
'Total' [1.8890e+003] [35] [] [] []
STATS =
source: 'anova2'
sigmasq: 50.8611
colmeans: [44.8889 40.6667 38 48.3333]
coln: 9
rowmeans: [42.2500 44.2500 42.4167]
rown: 12
inter: 1
pval: 0.9540
df: 24
同时返回图形,图形中显示了方差分析表。诸列(Columns)表示的是4种施肥方案,诸行(Rows)表示的是3种深翻方案,从上述方差分析表中可以看出,施肥方案有显著影响,深翻方案无显著影响,两因素间无交互作用。
COMPARISON = multcompare(STATS,'alpha',0.05)
结果显示
Note: Your model includes an interaction term. A test of main
effects can be difficult to interpret when the model includes
interactions.
COMPARISON =
1.0000 2.0000 -5.0520 4.2222 13.4964
1.0000 3.0000 -2.3853 6.8889 16.1631
1.0000 4.0000 -12.7187 -3.4444 5.8298
2.0000 3.0000 -6.6075 2.6667 11.9409
2.0000 4.0000 -16.9409 -7.6667 1.6075
3.0000 4.0000 -19.6075 -10.3333 -1.0591
多重比较的结果发现,第3,4种施肥方案差异显著。Matlab返回的注释说明,对于有交互项的情形,上述多重比较仅供参考。
5.2.2 无重复实验的双因素方差分析
要把交互作用与随机误差区别开,就必须对每种搭配进行重复试验。如果两个因素间确无交互作用,线性模型可以简化为:
,  (5-23)
此时检验H01H02即可。
利用anova2仍可进行方差分析,其原理仍是平方和分解与检验。不再一一罗列公式。
例5.4 某养猪场进行猪增重试验,选择4个品种的猪和3种饲料,共12中搭配方案,每种饲养一头,三个月后增重数据如表5-6所示。试研究品种与饲料对于猪增重的影响。
 复制粘贴输入数据矩阵x,执行
[P,TABLE,STATS] = anova2(x)
表5-6 猪增重数据表
 
饲料1
饲料2
饲料3
品种1
51
53
52
品种2
56
57
58
品种3
45
49
47
品种4
42
44
43
计算结果为:
P =
0.0156 0.0000
TABLE =
'Source' 'SS' 'df' 'MS' 'F' 'Prob>F'
'Columns' [ 10.5000] [ 2] [ 5.2500] [ 9] [ 0.0156]
'Rows' [332.2500] [ 3] [110.7500] [189.8571] [2.4683e-006]
'Error' [ 3.5000] [ 6] [ 0.5833] [] []
'Total' [346.2500] [11] [] [] []
STATS =
source: 'anova2'
sigmasq: 0.5833
colmeans: [48.5000 50.7500 50]
coln: 4
rowmeans: [52 57 47 43]
rown: 3
inter: 0
pval: NaN
df: 6
可以看出,诸列(Columns)间差异显著,说明饲料作用显著;诸行(Rows)间差异极显著,说明在这4个猪的品种间,增重差异极显著。
对于无重复实验可以可靠地进行多重比较。继续计算:
COMPARISON = multcompare(STATS,'alpha',0.05, 'estimate','column')
结果为
COMPARISON =
1.0000 2.0000 -3.9071 -2.2500 -0.5929
1.0000 3.0000 -3.1571 -1.5000 0.1571
2.0000 3.0000 -0.9071 0.7500 2.4071
发现第1,2种饲料差异显著。
COMPARISON = multcompare(STATS,'alpha',0.05, 'estimate','row')
结果为
COMPARISON =
1.0000 2.0000 -7.1588 -5.0000 -2.8412
1.0000 3.0000 2.8412 5.0000 7.1588
1.0000 4.0000 6.8412 9.0000 11.1588
2.0000 3.0000 7.8412 10.0000 12.1588
2.0000 4.0000 11.8412 14.0000 16.1588
3.0000 4.0000 1.8412 4.0000 6.1588
4个猪的品种,每种搭配对比,增重差异都显著。
COMPARISON = multcompare(STATS,'alpha',0.01, 'estimate','row')
结果为
COMPARISON =
1.0000 2.0000 -8.1014 -5.0000 -1.8986
1.0000 3.0000 1.8986 5.0000 8.1014
1.0000 4.0000 5.8986 9.0000 12.1014
2.0000 3.0000 6.8986 10.0000 13.1014
2.0000 4.0000 10.8986 14.0000 17.1014
3.0000 4.0000 0.8986 4.0000 7.1014
4个猪的品种,每种搭配对比,增重差异都极显著。
Read more ...