ff

1.变数与变量存档

1.1变量类型转换

1.1.1语法说明

$$double():转换为双精度$$

$$single():转换为单精度$$

$$int8():转换为8位有符号整数$$

$$int16():转换为16位有符号整数$$

$$uint8():转换为8位无符号整数$$

$$uint16():转换为16位无符号整数$$

1.1.2单一字元

定义：字符在ASCll中使用0到255之间的数字代码，表示通过将字符或字符串放入一对撇号中来创建字符或字符串。

ggghhhh

例：

s1 = 'h'
whos 
uint16(s1)

1.1.3多元字元

1.1.3.1说明

$$str=string(A):字符串数组$$

字符串数组而不是字符数组来表示文本。字符串数组的每个元素存储一个字符序列。序列可以具有不同长度，无需填充，例如 "yes" 和 "no"。只有一个元素的字符串数组也称为字符串标量。

按照标准数组运算对字符串数组进行索引、重构和串联，还可以使用 + 运算符向字符串追加文本。如果字符串数组表示数字，则可以使用 double 函数将其转换为数值数组。

1.1.3.2创建对象

创建字符串数组的一种方法是使用方括号将字符串串联成数组，就像将数字串联成数值数组一样。

例：

s1 = 'Example';
s2 = 'strring';
s3 = [s1 s2];
s4 = [s1;s2]

1.1.3.3逻辑运算和分配

许多数字和逻辑运算符可以应用于字符串

str = 'arrdvark'
'a' == str             

比较所有存在位置a的位置，存在出（1），不存在出（0）；若想改变a的存在，使a变为m即:

str(str == 'a') = 'm'

若想改变一串字符或者一句话的顺序即用：

s1='I like playing basketball'
s2=reverse(s1)

1.2结构体

1.2.1存储异构数据的方法

结构包含称为字段的数组，例：‘student’（结构）中分为四大类‘name’（姓名）、‘id’（身份件）、‘number’（数）、‘grade’（等级）。

将学生一的各项信息通过student进行数组存储。程序见下：

student.name = 'John Doe';
student.id = '500384197827647312';
student.number = 301073268;
student.grade = [100 75 73;95 91 85.5;100 98 72];
student

1.2.2向结构添加信息

存入信息不仅仅只限于一个数组的时候，往往都会添加另外一层结构，也就是‘student’分为多类。

添加学生二的各项信息程序只需要把‘student’表序号即可，具体程序如下：

 student(2).name = 'Ann Lane';
 student(2).id = '500384199836548712';
 student(2).number = 901673246;
 student(2).grade = [95 100 90;95 91 87;97 99 82];
 student

2.2.2.1延伸拓展

将已经存储好的2位同学，如果我想要知道第二位同学的

grade中的第3行第二列成绩，程序如下：

student(2).grade(3,2)

1.3结构功能

1.3.1功能基础语法

$$structArray=cell2struct(cellArray,fields,dim)$$

将元胞数组转换为结构体数组.通过元胞数组 cellArray 中包含的信息创建一个结构体数组 structArray。

fields 参数指定结构体数组的字段名称。

dim 参数向 MATLAB 指示创建结构体数组时要使用的元胞数组的轴。
$$fields=fieldnames(S)$$
返回元胞数组中结构体数组 S 的字段名称。
$$s=rmfield(s,field)$$
从结构体数组 s 中删除指定的一个或多个字段。使用字符向量元胞数组或字符串数组指定多个字段。s 的维度保持不变。

其中field为字段名称，指定为字符数组、字符向量元胞数组或字符串数组。

1.3.2功能语法举例

根据**1.2.2向结构添加信息**所构造出的例子student已记录了两名学生的各项信息，现将student中id信息抓取出来并且删除，操作视图如下：

fieldnames(student)
rmfield(student,'id')

1.4嵌套结构

以下是简单的嵌套函数结构

其运行代码为：

A=struct('data',[3 4 7;8 0 1],'nest', ...
    struct('testman','Test 1', ...
    'stare',[4 2 6],'ydata',[7 1 6]));
A(2).data=[9 2 3; 7 6 5];
A(2).nest.testnum='Test 2';
A(2).nest.xdata=[3 4 2];
A(2).nest.ydata=[5 0 9];
A.nest

1.5元胞数组

元胞数组存储异构数据的另一种方法，类似于矩阵，但每个条目包含不同类型的数据，使用{}声明。例如：

A(1,1)={[1 4 3;0 5 8;7 2 9]};
A(1,2)={'Anne Smith'};
A(2,1)={3+7i};
A(2,2)={-pi:pi:pi};
A

或者

A{1,1}=[1 4 3;0 5 8;7 2 9];
A{1,2}='Anne Smith';
A{2,1}=3+7i;
A{2,2}=-pi:pi:pi;
A

1.5.1访问单元阵列

大括号{}用于访问单元格arry的“内容”，要精确到一个具体的位置就是用矩阵来表示出位置，位于几行几列(,)例如：

 A{1,1}=[1 4 3;0 5 8;7 2 9];
A{1,2}='Anne Smith';
A{2,1}=3+7i;
A{2,2}=-pi:pi:pi;
A{1,1}(3,1)       %找出我需要的[3;1]数字

1.6单元格数组函数

1.6.1语法基本说明

$$C=cell(n)返回由空矩阵构成的n\times n元胞数组。$$

元胞数组是一种包含名为元胞的索引数据容器的数据类型，其中的每个元胞都可以包含任意类型的数据。元胞数组通常包含文本列表、文本和数字的组合或者不同大小的数值数组。通过将索引括在圆括号 () 中可以引用元胞集。使用花括号 {} 进行索引来访问元胞的内容。
$$C=num2cell(A)$$
通过将 A 的每个元素放置于 C 的一个单独元胞中，来将数组 A 转换为元胞数组 C。
$$C=mat2cell(A,rowDist)$$
将数组 A 划分为一个 n×1 元胞数组 C，其中 n 等于 rowDist 中元素的数量。

1.6.2数组函数举例

把amatrix变换成一个单元格变量，把结构体里的数据变成cell的数组，具体流程如下：

实现关键代码num2cell，例：

a=magic(3);
b=num2cell(a);
c=mat2cell(a,[1 1 1],3)

1.7多维数组

1.7.1多维数组第一种编辑方式

构造一个3维数组，将每一个数组标记位置，宣告方式{row,column,layer}，例：

%3维数组
A{1,1,1}=[1 2;4 5];
A{1,2,1}='Name';
A{2,1,1}=2-4i;
A{2,1,2}=7;
A{1,1,2}='Name2';
A{1,2,2}=3;
A{2,1,2}=0:1:3;
A{2,2,2}=[4 5]';
A

1.7.2多维数组第二种编辑方式

1.7.2.1串联数组语法

$$C=cat(dim,A1,A2,\dots ,An)$$

沿维度 dim 串联 A1、A2、…、An。

可以使用方括号运算符 [] 进行串联。例如，[A,B] 或 [A B] 将水平串联数组 A 和 B，而 [A; B] 将垂直串联它们。

1.7.2.2使用示例

数组串联，类似于高斯消去法用法，运作效果如下：

A=[1 2;3 4];
B=[5 6;7 8];
C=cat(1,A,B);
C=cat(2,A,B);
C=cat(3,A,B);
C

或许先建成一个二维的各项数组，再用cat进行一个数组统合，形成大的结构性数组，例：

A{1,1} = [1 2;4 5];
A{1,2} = 'Name';
A{2,1} = 2-4i;
A{2,2} = 7;
B{1,1} = 'NAME2';
B{1,2} = 3;
B{2,1} = 0:1:3;
b{2,2} = [4,5]';
C = cat(3,A,B)

1.8重构数组

1.8.1重构基础语法

$$B=reshape(A,sz1,...,szN)$$

将 A 重构为一个 sz1×...×szN 数组，其中 sz1,...,szN 指示每个维度的大小。可以指定 [] 的单个维度大小，以便自动计算维度大小，以使 B 中的元素数与 A 中的元素数相匹配。

1.8.2使用方式

返回具有指定的rouw和列的nw数组，例：

A = {'James Bond', [1 2;3 4;5 6]; pi, magic(5)}
C = rehsape(A,1,4)

C是对A的数组还原，以现有A的数据构建一个关于C的单维数组

1.9检查bariable和variable状态

检查一个数组是否具有cell、real、char等情况，需要用到的命令：
$$isinteger:确定输入是否为整数数组$$

$$isreal:确定数组是否为实数数组$$

$$iscell:确定输入是否为元胞数组$$

$$ischar:确定输入是否为字符数组$$

$$isempty:确定数组是否为空$$

以上都是常见的使用状态，具体使用可以在matlab中使用help (函数命令)进行查询使用。

2.文件存取

2.1文件的高阶存储方法

2.1.1save()和load()

将（全部）wrokspace数据保存到文件,附加添加功能，例：

a=magic(4);
save mydata1.mat
save mydata2.mat -ascii

加载存储在文件中的数据：

load('mydata.mat')
load('mydata2.mat','-ascii')

2.2 Excel信息导入即文件修订

2.2.1导入语法

$$num=xlsread(filename,sheet,{}^{\prime }{m}^{\prime },1,xlRange)$$

读取指定的工作表和范围。·m·更改或者添加的数值，1表示m添加的位置xlRange。

2.2.2 xlsread()举例

data=xlsread('D:/高等数学代码/原始数据.xlsx');
M=mean(data')';
xlswrite('原始数据xlsx',M,1,'F2:F13')

此代码将已知数据算出平均值，然后加到’F2:F13’中。

2.3低阶存储方法

2.3.1低级文件l/o函数说明

$$fopen:打开文件或获得有关打开文件的信息$$

$$fopen(filename{,}^{\prime }permissio{n}^{\prime })$$

打开文件 filename 以便以二进制读取形式进行访问，并返回等于或大于 3 的整数文件标识符。

permission常用的几个读取方式

‘r’	‘r+’	‘w’
’w+'	’a’	’a+'

最常用的就是r和w。
$$fprintf(fileID,formatSpec,A1,...,An)$$
按列顺序将 formatSpec 应用于数组 A1,...An 的所有元素，并将数据写入到一个文本文件。fprintf 使用在对 fopen 的调用中指定的编码方案。

2.3.2 l/o打开和关闭文件

x=0:pi/10:pi;y=sin(x);         %生成x,y
fid=fopen('sinx.txt','w');     %打开文件
for i=1:11
    fprintf(fid,'%5.3f %8.4f\n',x(i),y(i));  %把x,y写入文件
end
fclose(fid);                   %关闭文件
type sinx.txt                  %打开文件内容

再去对应运行的文件夹下面找到’sinx.txt’文件，打开就可以看见其中的内容。

2.3.3通过格式化的l/o读写

2.3.3.1读取函数fscanf()

$$A=fscanf(fileID,formatSpec,sizeA)$$

将文件数据读取到维度为 sizeA 的数组 A中，并将文件指针定位到最后读取的值之后。fscanf 按列顺序填充 A。

formatSpec 指定的格式解释文件中的值。

sizeA 必须为正整数或采用 [m n] 的形式，其中 m 和 n 为正整数。

2.3.3.2 示例读取

根据txt所给的信息如下，用matlab进行读取。

John	1995	12	5	12.3	3.24
Tom	1995	12	7	2.3	2
Jean	1996	3	2	10.2	0

运行程序如下：

clc;clear
fid = fopen('files.txt','r');
i=1;
while ~feof(fid)
    name(i,:)=fscanf(fid,'%5c',1);
    year{i}=fscanf(fid,'%d',1);
    no1{i}=fscanf(fid,'%d',1);
    no2{i}=fscanf(fid,'%d',1);
    no3{i}=fscanf(fid,'%g',1);
    no4{i}=fscanf(fid,'%g\n',1);
    i=i+1;
end
fclose(fid)

2.4拓展延伸数值字段

数值字段类型	转换设定符	详细信息
有符号整数	%d	以 10 为基数
%i	文件中的值确定相应基数：默认值以 10 为基数。如果初始数字为 0x 或 0X，则值为十六进制（以 16 为基数）。如果初始数字为 0，则值为八进制（以 8 为基数）。
%ld 或 %li	64 位值，以 10、8 或 16 为基数
无符号整数	%u	以 10 为基数
%o	以 8 为基数（八进制）
%x	以 16 为基数（十六进制）
%lu、%lo、%lx	64 位值，以 10、8 或 16 为基数
浮点数	%f	浮点字段可以包含下列任意项（不区分大小写）：Inf、-Inf、NaN 或 -NaN。

3.初阶绘图

3.1图线的绘制与装饰

3.1.1 plot函数绘制图线

$$plot(X1,Y1,LineSpec1,...,Xn,Yn,LineSpecn)$$

设置每个线条的线型、标记符号和颜色。

• 如果 X 和 Y 都是向量，则它们的长度必须相同。plot 函数绘制 Y 对 X 的图。
• 如果 X 和 Y 均为矩阵，则它们的大小必须相同。plot 函数绘制 Y 的列对 X 的列的图。
• 如果 X 或 Y 中的一个是向量而另一个是矩阵，则矩阵的各维中必须有一维与向量的长度相等。如果矩阵的行数等于向量长度，则 plot 函数绘制矩阵中的每一列对向量的图。如果矩阵的列数等于向量长度，则该函数绘制矩阵中的每一行对向量的图。如果矩阵为方阵，则该函数绘制每一列对向量的图

线型符号	线型设定符	标记	标记设定符	颜色	颜色设定符
-	实线（默认）	o	圆圈	y	黄色
--	虚线	+	加号	m	品红色
:	点线	*	星号	c	青蓝色
-.	点划线	.	点	r	红色
		x	叉号	g	绿色
		s	方形	b	蓝色
		d	菱形	w	白色
		^	上三角	k	黑色
		v	下三角
		>	右三角
		<	左三角
		p	五角形
		h	六角形

例如：绘制一个(0,2π)内余弦函数的图像:

x = 0:pi/20:2*pi;
y = cos(x);
plot(x, y, 'r.-')

3.1.2装饰图线

使用legend()函数为图片增加图例

使用title()和*label()为图片增加标题和标签

x = 0:0.1:2*pi; y1 = sin(x); y2 = exp(-x);
plot(x, y1, '--*', x, y2, ':o');
xlabel('t = 0 to 2\pi');
ylabel('values of sin(t) and e^{-x}')
title('Function Plots of sin(t) and e^{-x}');
legend('sin(t)','e^{-x}');

使用text()和annotation()为图片增加注解

x = linspace(0,3); y = x.^2.*sin(x); plot(x,y);
line([2,2],[0,2^2*sin(2)]);
str = '$$ \int_{0}^{2} x^2\sin(x) dx $$';
text(0.25,2.5,str,'Interpreter','latex');
annotation('arrow','X',[0.32,0.5],'Y',[0.6,0.4]); 

3.2控制坐标轴,边框与网格

$$subplot(m,n,p)$$
将当前图窗划分为 m×n 网格，并在 p 指定的位置创建坐标区。MATLAB 按行号对子图位置进行编号。第一个子图是第一行的第一列，第二个子图是第一行的第二列，依此类推。如果指定的位置已存在坐标区，则此命令会将该坐标区设为当前坐标区。

使用下列命令可以控制坐标轴,边框与网格.

命令	作用
grid on/off	设置网格可见性
box on/off	设置边框可见性
axis on/off	设置坐标轴可见性
axis normal	还原默认行为,将图框纵横比模式和数据纵横比模式的属性设置为自动
axis square	使用相同长度的坐标轴线,相应调整数据单位之间的增量
axis equal	沿每个坐标轴使用相同的数据单位长度
axis tight	将坐标轴范围设置为等同于数据范围,使轴框紧密围绕数据

例：演示改变圆形的X,Y的4类不同情况。

t = 0:0.1:2*pi; x = 3*cos(t); y = sin(t);
subplot(2, 2, 1); plot(x, y); axis normal
subplot(2, 2, 2); plot(x, y); axis square
subplot(2, 2, 3); plot(x, y); axis equal
subplot(2, 2, 4); plot(x, y); axis equal tight

3.3绘制多条图线

3.3.1在一个图像上绘制多条图线

默认情况下,每次执行plot()函数都会清除上一次绘图的结果,多次执行plot()只会保留最后一次绘制的图形.

plot(cos(0:pi/20:2*pi));
plot(sin(0:pi/20:2*pi));

我们可以使用hold on和hold off命令控制绘图区域的刷新,使得多个绘图结果同时保留在绘图区域中.

hold on		% 提起画笔,开始绘制一组图片			
plot(cos(0:pi/20:2*pi));
plot(sin(0:pi/20:2*pi));
hold off	% 放下画笔,该组图片绘制完毕

3.3.2在一个窗口内绘制多个图像

使用subplot()函数可以在一个窗口内绘制多个图像.

subplot(2,2,1);
x = linspace(-3.8,3.8);
y_cos = cos(x);
plot(x,y_cos);
title('Subplot 1: Cosine')

subplot(2,2,2);
y_poly = 1 - x.^2./2 + x.^4./24;
plot(x,y_poly,'g');
title('Subplot 2: Polynomial')

subplot(2,2,[3,4]);
plot(x,y_cos,'b',x,y_poly,'g');
title('Subplot 3 and 4: Both')

3.4图形对象的操作

3.4.1获取图形句柄

图形句柄本质上就是一个浮点数,可以唯一确定一个图形对象.下面几个函数用于获取图形句柄.

Function	Purpose
gca()	获取当前坐标轴的句柄
gcf()	获取当前图像的句柄
allchild(handle_list)	获取该对象的所有子对象的句柄
ancestor(h,type)	获取对象最近的type类型的祖先节点
delete(h)	删除某对象
findall(handle_list)	获取该对象的后代对象

所有绘图函数也会返回图形对象的句柄.

3.4.2通过图形句柄操作图形属性

使用get()和set()函数可以对图形对象的属性进行访问和修改.
$$v=get(h,propertyName)$$
返回特定属性 propertyName 的值。使用时须用单引号将属性名引起来，例如，get(h,'Color')。如果您不指定输出参数，则 MATLAB 会在屏幕上显示该信息。
$$set(H,Name,Value)$$
为 H 标识的对象指定其 Name 属性的值。使用时须用单引号将属性名引起来，例如，set(H,'Color','red')。如果 H 是对象的向量，则 set 会为所有对象设置属性。

3.5将图形保存到文件

$$saveas(fig,filename)$$

使用saveas(fig,filename)命令可以将图形对象保存到文件中,其中fig为图形句柄,filname为文件名.

例如 'myplot.jpg'。文件扩展名用于定义文件格式。如果不指定扩展名，则 saveas 会将图窗保存为 FIG 文件。要保存当前图窗，请将 fig 指定为 gcf

使用saveas()函数将图像保存成位图时,会发生失真.要精确控制生成图片的质量,可以使用print()函数,
$$print(filename,formattype)$$
使用指定的文件格式将当前图窗保存到文件中，例如 print('BarPlot','-dpng')。如果该文件不包括扩展名，则 print 会附加适用的扩展名。

4.进阶绘图

4.1二维图表

4.1.1折线图

函数	图形描述
loglog()	x轴和y轴都取对数坐标
semilogx()	x轴取对数坐标,y轴取线性坐标
semilogy()	x轴取线性坐标,y轴取对数坐标
plotyy()	带有两套y坐标轴的线性坐标系
ploar()	极坐标系

4.1.1.1对数坐标系

演示对数坐标系图线：

clc;clear
x = logspace(-1,1,100);       
y = x.^2;
subplot(2,2,1);
plot(x,y)
title('plot');
subplot(2,2,2);
semilogx(x,y);
title('semilogx');
subplot(2,2,3);
semilogy(x,y);
title('semilogy');
subplot(2,2,4)
loglog(x,y);
title('loglog')

演示结果：

4.1.2双y轴图线

$$plotyy(X1,Y1,X2,Y2)$$

制 Y1 对 X1 的图，在左侧显示 y 轴标签，并同时绘制 Y2 对 X2 的图，在右侧显示 y 轴标签。

演示双y轴图线：

x = 0:0.01:20;
y1 = 200*exp(-0.05*x).*sin(x);
y2 = 0.8*exp(-0.5*x).*sin(10*x);

figure % new figure
[hAx,hLine1,hLine2] = plotyy(x,y1,x,y2);

title('Multiple Decay Rates')
xlabel('Time (\musec)')

ylabel(hAx(1),'Slow Decay') % 左边 y-axis 
ylabel(hAx(2),'Fast Decay') % 右边 y-axis

演示结果如下：

4.1.3统计图表

4.3.1正态分布确认随机标量

$$X=randn返回一个从标准正态分布中得到的随机标量。$$

$$X=randn(n)返回由正态分布的随机数组成的n\times n矩阵。$$

4.1.3.2直方图

$$hist(x)$$

基于向量 x 中的元素创建直方图条形图。x 中的元素有序划分入 x 轴上介于 x 的最小值和最大值间的 10 个等间距 bin 中。hist 将 bin 显示为矩形，这样每个矩形的高度就表示 bin 中的元素数量。

如果输入是多列数组，则 hist 为每列 x 创建直方图并将它们叠加到一个绘图上。

如果输入为 categorical 数据类型，则每个 bin 是一个 x 类别。

演示直方图：

clc;clear;close all
y = randn(1,1000);
subplot(2,1,1);
hist(y,10);
title('Bins = 10');
subplot(2,1,2);
hist(y,50);
title('Bins = 50');

演示结果如下：

4.4条形图

$$bar(y)$$

创建一个条形图，y 中的每个元素对应一个条形。如果 y 是矩阵，则 bar 根据 y 中的行对条形分组。
$$bar3(Z)$$
绘制三维条形图，Z 中的每个元素对应一个条形图。如果 Z 是向量，y 轴的刻度范围是从 1 至 length(Z)。如果 Z 是矩阵，则 y 轴的刻度范围是从 1 到 Z 的行数。

演示条形图1：

clc;clear;close all
x = [1 2 5 4 8];   y = [x;1:5];
subplot(1,3,1); bar(x); title('A bargraph of vector x')
subplot(1,3,2); bar(y); title('A bargraph of vector y')
subplot(1,3,3); bar3(y); title('A 3D bargraph')

演示结果1：

演示条形图2：

clc;clear;close all
x = [1 2 5 4 8];   y = [x;1:5];
subplot(1,2,1);
bar(y,'stacked');
title('stacked');
subplot(1,2,2);
barh(y,'stacked');
title('Horizontal')

演示结果2：

4.5饼图

$$pie(X)$$

使用 X 中的数据绘制饼图。饼图的每个扇区代表 X 中的一个元素。

• 如果 sum(X) ≤ 1，X 中的值直接指定饼图扇区的面积。如果 sum(X) < 1，pie 仅绘制部分饼图。

• 如果 sum(X) > 1，则 pie 通过 X/sum(X) 对值进行归一化，以确定饼图的每个扇区的面积。

• 如果 X 为 categorical 数据类型，则扇区对应于类别。每个扇区的面积是类别中的元素数除以 X 中的元素数的结果。
  $$pie3(X)$$
  使用 X 中的数据绘制三维饼图。X 中的每个元素表示饼图中的一个扇区。

演示饼图：

clc;clear；close all
a = [10 5 20 30];
subplot(1,3,1);  pie(a);
subplot(1,3,2);  pie(a,[1,1,1,1]);
subplot(1,3,3);  pie3(a,[0,0,0,1]);

4.6极坐标图线

$$polar(theta,rho)$$

创建角 theta 对半径 rho 的极坐标图。theta 是从 x 轴到半径向量所夹的角（以弧度单位指定）；rho 是半径向量的长度（以数据空间单位指定）。

演示极坐标图形：

clc;clear;close all
x = 1:100; theta = x/10;  r = log10(x);
subplot(2,2,1);  
polar(theta,r); %螺旋线
theta = linspace(0,2*pi); 
r = cos(4*theta);
subplot(2,2,2); 
polar(theta,r);%花瓣
theta = linspace(0,2*pi,6);
r = ones(1,length(theta));
subplot(2,2,3); 
polar(theta,r);%五边形
theta = linspace(0,2*pi);
r = 1-sin(theta);
subplot(2,2,4); 
polar(theta,r);%心形线

演示结果：

4.7阶梯图和离散序列数据

$$stairs(Y)$$

绘制 Y 中元素的阶梯图。

• 如果 Y 为向量，则 stairs 绘制一个线条。
• 如果 Y 为矩阵，则 stairs 为每个矩阵列绘制一个线条。

$$stem(Y)$$

将数据序列 Y 绘制为从沿 x 轴的基线延伸的针状图。各个数据值由终止每个针状图的圆指示。

• 如果 Y 是向量，x 轴的刻度范围是从 1 至 length(Y)。
• 如果 Y 是矩阵，则 stem 将根据相同的 x 值绘制行中的所有元素，并且 x 轴的刻度范围是从 1 至 Y 中的行数。

演示阶梯图和离散序列数据：

clc;clear;close all
x=linspace(0,4*pi,40);
y=sin(x);
subplot(1,2,1);
stairs(y);
subplot(1,2,2);
stem(y)

演示结果：

4.8方框图

$$boxplot(x)$$

x）创建x中数据的方框图。如果x是矢量，则方框图绘制一个方框。如果x是一个矩阵，那么boxplot为x的每一列绘制一个框。

演示方框图：

clc;clear;close all
load carsmall
boxplot(MPG,Origin)

演示结果：

4.9含误差条的线图

$$errorbar(y,err)$$

创建 y 中数据的线图，并在每个数据点处绘制一个垂直误差条。err 中的值确定数据点上方和下方的每个误差条的长度，因此，总误差条长度是 err 值的两倍。

演示含误差条的线图：

clc;clear;close all
x=0:pi/10:pi;
y=sin(x);
e=std(y)*ones(size(x));
errorbar(x,y,e);

演示结果：

4.2三维图表

4.2.1二维图与三维图间的关系

4.2.1.1二维图转为三维图

在MATLAB中,所有的图都是三维图,二维图只不过是三维图的一个投影.点击图形窗口的Rotate 3D按钮,即可通过鼠标拖拽查看该图形的三维视图.

4.2.1.2三维图转为二维图

$$imagesc(C)$$

将数组 C 中的数据显示为一个图像，该图像使用颜色图中的全部颜色。C 的每个元素指定图像的一个像素的颜色。生成的图像是一个 m×n 像素网格，其中 m 和 n 分别是 C 中的行数和列数。这些元素的行索引和列索引确定了对应像素的中心。

演示三维图转为二维图：

clc;clear;close all
[x,y]=meshgrid(-3:.2:3,-3:.2:3);
z=x.^2+x.*y+y.^2;
surf(x,y,z);
box on
set(gca,'FontSize',16);
zlabel('x');
xlim([-4 4]);  xlabel('x');  ylim([-4 4]);  ylabel('y')
figure(1)
help on
figure(2)
imagesc(z);   axis square;
xlabel('x');   ylabel('y')

演示结果：

三维转变二维

补充：使用colorbar命令可以在生成的二维图上增加颜色与高度间对应关系的图例,使用colormap命令可以改变配色方案.

4.2.2三维图的绘制

4.2.2.1生成二维网格

$$meshgrid(x,y)$$

基于向量 x 和 y 中包含的坐标返回二维网格坐标。X 是一个矩阵，每一行是 x 的一个副本；Y 也是一个矩阵，每一列是 y 的一个副本。坐标 X 和 Y 表示的网格有 length(y) 个行和 length(x) 个列。

meshgrid()函数将输入的两个向量进行相应的行扩充和列扩充以得到两个增广矩阵,对该矩阵可应用二元函数.

演示生成二维网格矩阵：

clc;clear
x=-2:1:2;
y=-2:1:2;
[x,y]=meshgrid(x,y);
z=x.^2+y.^2

演示结果：

x =

-2    -1     0     1     2
-2    -1     0     1     2
-2    -1     0     1     2
-2    -1     0     1     2
-2    -1     0     1     2

y =

-2    -2    -2    -2    -2
-1    -1    -1    -1    -1
 0     0     0     0     0
 1     1     1     1     1
 2     2     2     2     2

z =

 8     5     4     5     8
 5     2     1     2     5
 4     1     0     1     4
 5     2     1     2     5
 8     5     4     5     8

4.2.2.2绘制三维线

$$plot3(X,Y,Z)$$

绘制三维空间中的坐标。

• 要绘制由线段连接的一组坐标，请将 X、Y、Z 指定为相同长度的向量。
• 要在同一组坐标轴上绘制多组坐标，请将 X、Y 或 Z 中的至少一个指定为矩阵，其他指定为向量。

演示三维线：

clc;clear;close all
x=0:0.1:3*pi;  z1=sin(x);   z2=sin(2.*x);  z3=sin(3.*x);
y1=zeros(size(x));  y3=ones(size(x));   y2=y3./2;
plot3(x,y1,z1,x,y2,z2,x,y3,z3);
grid on
xlabel('x-axis');    ylabel('y-axis');    zlabel('z-axis');

演示结果：

4.2.2.3绘制三维表面图

绘制三维表面图有两种方法，一种使用mesh另外一种使用surf。mesh不会填充网格，surf会填充网格。
$$mesh(X,Y,Z)$$
使用 Z 确定的颜色绘制线框网格，因此其颜色与曲面高度成比例。如果 X 和 Y 为向量，length(X) = n 且 length(Y) = m，其中 [m,n] = size(Z)。在本示例中，(X(j), Y(i), Z(i,j)) 是线框网格线的交点；X 和 Y 分别对应于 Z 的列和行。如果 X 和 Y 为矩阵，则 (X(i,j), Y(i,j), Z(i,j)) 是线框网格线的交点。X、Y 或 Z 中的值可以是数值、日期时间值、持续时间值或分类值。
$$surf(X,Y,Z)$$
创建一个三维曲面图。该函数将矩阵 Z 中的值绘制为由 X 和 Y 定义的 x-y 平面中的网格上方的高度。函数还对颜色数据使用 Z，因此颜色与高度成比例。

演示绘制三维表面图：

clc;clear;close all
x=-3.5:0.2:3.5;
y=-3.5:0.2:3.5;
[x,y]=meshgrid(x,y);
z=x.*exp(-x.^2-y.^2);
subplot(1,2,1);   mesh(x,y,z)
subplot(1,2,2);   surf(x,y,z)

演示结果：

4.2.2.4绘制三维图形的等高线

绘制三维表面图有两种方法，一种使用contour另外一种使用contourf。contour不会填充网格，contourf会填充网格。
$$contour(Z)$$
创建一个包含矩阵 Z 的等值线的等高线图，其中 Z 包含 x-y 平面上的高度值。MATLAB® 会自动选择要显示的等高线。Z 的行索引和列索引分别是平面中的 x 和 y 坐标。
$$contourf(Z)$$
创建一个包含矩阵 Z 的等值线的填充等高线图，其中 Z 包含 x-y 平面上的高度值。MATLAB® 会自动选择要显示的等高线。Z 的行索引和列索引分别是平面中的 x 和 y 坐标。

演示绘制三维图形的等高线：

clc;clear;close all
x=-3.5:0.2:3.5;
y=-3.5:0.2:3.5;
[x,y]=meshgrid(x,y);
z=x.*exp(-x.^2-y.^2);
subplot(2,2,1:2);   
mesh(x,y,z);
axis square;
subplot(2,2,3); 
contour(x,y,z);
axis square;    
subplot(2,2,4);  
contourf(x,y,z);
axis square;

演示结果：

4.2.2.5三维图的视角与打光

4.2.2.5.1调整角度

$$view(v)$$

根据 v（二元素或三元素数组）设置视线：

• 二元素数组 - 其值分别是方位角和仰角。
• 三元素数组 - 其值是从图框中心点到照相机位置所形成向量的 x、y 和 z 坐标。MATLAB® 使用指向同一方向的单位向量计算方位角和仰角。

view()函数接受两个浮点型参数,分别表示两个方位角azimuth和elevation.

演示调整角度：

clc;clear;close all
sphere(50);  shading flat;
light('position',[1 3 2]);
light('position',[-3 -1 3]);
material shiny;
axis vis3d off;
set(gcf,'color',[1 1 1])
view(-45,20)

演示结果：

4.2.2.5.2调整灯光

$$light{(}^{\prime }PropertyNam{e}^{\prime },propertyvalue,...)$$

使用给定属性的指定值创建一个 Light 对象。有关上述属性的说明，请参阅 Light 属性。MATLAB®软件将该光源指定为当前坐标区的父级，除非您使用 Parent 属性指定了其他坐标区。

画出红色圆球：

clc;clear;close all
[x,y,z]=sphere(64); h=surf(x,y,z);
axis square vis3d off;
reds=zeros(256,3); reds(:,1)=(0:256.-1)/255;
colormap(reds); shading interp; lighting phon

演示调整灯光:

set(h,'AmbientStrength',0.75,'DiffuseStrength',0.5);
L1=light('position',[-1,-1,-1]);
set(L1,'position',[-1,-1,-1]);

演示结果：

4.2.2.6绘制三维体

$$patch(X,Y,C)$$

使用 X 和 Y 的元素作为每个顶点的坐标，以创建一个或多个填充多边形。patch 以您指定顶点的顺序连接这些顶点。要创建一个多边形，请将 X 和 Y 指定为向量。要创建多个多边形，请将 X 和 Y 指定为矩阵，其中每一列对应于一个多边形。C 决定多边形的颜色。

演示绘制三维体：

v=[0 0 0;1 0 0;1 1 0;0 1 0;0.25 0.25 1;...
    0.75 0.25 1;0.75 0.75 1;0.25 0.75 1];
f=[1 2 3 4;5 6 7 8;1 2 6 5;2 3 7 6;3 4 8 7;4 1 5 8];
subplot(2,2,[1,3]);patch('Vertices',v,'Faces',f,...
    'FaceVertexCData',hsv(6),'FaceColor','flat');
view(3);axis square tight; grid on;
subplot(2,2,[2,4]);patch('Vertices',v,'Faces',f,...
    'FaceVertexCData',hsv(8),'FaceColor','interp');
view(3);axis square tight;  grid on

演示结果：

5.图形界面—GUI设计

5.1打开设计GUI

在MATLAB行窗命令口中输入guide会出现以下界面：

点击Blank Gul (Ddfault)会出现以下界面：

5.2设计第一个简单的程序

在untitled1.fig中选取3个普通按钮以及两个坐标轴，并使用工具中的对齐现象进行对齐,以及利用Fontisize调整字体和String改变标题 ，保存。例：

打开MATLAB中untitled1.m找到fuction中的pushbutton1..Callback在最后一部分加入我们想表达的图像，比如：

handles.peaks=peaks(35);
handles.membrane=membrane;
[x,y]=meshgrid(-8:.5:8);
r=sqrt(x.^2+y.^2)+eps;
sinc=sin(r)./r;
handles.sinc=sinc;
handles.current_data=handles.peaks;
axes(handles.axes1);
surf(handles.current_data)

保存，并运行，形成如下情况：

5.3设计第二个简单程序

设计一个A+B=？，其中为A一个单独的一个未知数范围（0—100），B同为一个未知数范围（0——100）。

创建一个GUI并命名**（untitled，fig），准备五个静态文本**，利用Fontisize调整字体和String改变标题 ，并排版好，插入2个滑动条位于2个静态文本的中心位置，改变范围（0—100）。例：

打开MATLAB中untitled.m找到fuction中的slider2_Callback在最后一部分加入A与B的关键定义，以及逻辑运算，比如：

A = get(handles.slider1,'Value');
B = get(handles.slider2,'Value');
SUM = int16(A+B);
set(handles.text2,'String',['A+B=',num2str(SUM)])

保存，并运行，形成如下情况：

5.4 存储变成档案

$$deploytool$$

deploytool打开编译器应用程序的列表。找到指定的GUI,把GUI进行打包，compile（链接），以程序的方式进行储存。

可以把这段程序带到其他电脑进行运行。

6.影像处理

6.1图像在MATLAB中的储存格式

MATLAB能够处理的数字图像分为三种:二值图像,灰度图像,彩色图像.

二值图像在MATLAB中以一个矩阵存储,矩阵中元素的取值为0(表示白)或1(表示黑).

灰度图像在MATLAB中以一个矩阵存储,矩阵中元素的取值介于0~255之间,表示灰度.

彩色图像在MATLAB中以三个矩阵存储,每个矩阵中元素的取值介于0~255之间,分别表示颜色R,G,B分量的浓度

6.2读取和展示图像

$$A=imread(filename)$$

从 filename 指定的文件读取图像，并从文件内容推断出其格式。如果 filename 为多图像文件，则 imread 读取该文件中的第一个图像。

使用imread()函数将图像读取到内存中,使用imshow()函数展示图像,使用imwrite()函数将内存中的图像写进硬盘.

演示读取图像：

clc;clear;close all
I=imread('pout.tif');%read
 imshow(I);%show

演示结果：

补充：可以使用imageinfo继续查看图片文件的详细信息；可以使用imtool可以打开图像处理工具。

6.3图像的运算

6.3.1图像的四则运算

要对两个图像进行四则运算,要求这两个图像的尺寸相同.下面是常用的图像四则运算函数.

函数	作用
imabsdiff()	两个图像求差值
imadd()	一个图像加上另一个图像或常数
imsubtract()	一个图像减去另一个图像或常数
immultiply()	一个图像乘以另一个图像或常数
imdivide()	一个图像除以另一个图像或常数
imcomplement()	对图像取反

例：演示两个图像差值情况

clc;clear;close all
I=imread('rice.png');
subplot(2,2,[1,3]);   imshow(I);
J=immultiply(I,1.5);
subplot(2,2,[2,4]);   imshow(J)

演示结果：

[注：左图像的各个像素乘以1.5倍变成右图像]

例：演示两个图像相加的情况

clc;clear;close all
I=imread('rice.png');
J=imread('cameraman.tif');   K=imadd(I,J);
subplot(2,3,[1,4]);    imshow(I);
subplot(2,3,[2,5]);    imshow(K);
subplot(2,3,[3,6]);    imshow(J);

演示结果：

[注:可以看到,进行加法操作后,得到的图像比原本的两个都亮,这是因为图像矩阵的数值整体上增加了.]

6.3.2像素的统计分布

$$[counts,binLocations]=imhist(I)$$

计算灰度图像I的直方图。imhist函数返回以计数为单位的直方图计数和以二进制位置为单位的二进制位置。直方图中的箱数由图像类型确定。

演示灰度图像：

clc;clear;close all
I=imread('pout.tif');
imhist(I)

演示结果：

6.3.2.1图像均衡化

$$J=histeq(I,hgram)$$

灰度图像I使得具有长度（hgram）bin的输出灰度图像J的直方图与目标直方图hgram近似匹配。

演示图像均衡化前与均衡化后对比：

clc;clear;close all
I=imread('pout.tif');  I2=histeq(I);
subplot(2,4,[1,5]);   imhist(I);
subplot(2,4,[2,6]);   imshow(I);
subplot(2,4,[3,7]);   imshow(I2);
subplot(2,4,[4,8]);   imhist(I2);

演示结果：

6.3.2.2图像的二值化

将灰度图像变为二值图像的过程被称为二值化,MATLAB内置了两个与二值化相关的函数.
$$T=graythresh(I)$$
选择一个阈值，使阈值化的黑白像素的组内方差最小化。全局阈值T可与imbinarize一起用于将灰度图像转换为二值图像。

利用in2bw进行二值化变化.

clc;clear;close all
I = imread('rice.png'); 
level=graythresh(I); bw=im2bw(I, level); 
figure(1);imhist(I); 
figure(2);
subplot(1,2,1); imshow(I);
subplot (1,2,2); imshow(bw)

演示结果：

6.4图像的几何变换

图像的几何变换本质上就是将图像乘以一个矩阵得到新图像的过程.

变换形式	图形示意	数学变换	MATLA命令
位移(Translation)		$$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ 1\end{array}\right]=\left[\begin{array}{ccc}1& 0& t_x\\ 0& 1& t_y\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]$$	imtranslate()
缩放(Scale)		$$\left[\begin{array}{ccc}{s}_x& 0& 0\\ 0& s_x& 0\\ 0& 0& 1\end{array}\right]$$	imresize()
错切(Shear)		$$\left[\begin{array}{ccc}1& h_x& 0\\ h_y& 1& 0\\ 0& 0& 1\end{array}\right]$$
旋转(Rotate)		$$\left[\begin{array}{ccc}cos\theta & sin\theta & 0\\ -sin\theta & cos\theta & 0\\ 0& 0& 1\end{array}\right]$$	imrotate()

演示图像变化：

clc;clear;close all
I=imread('rice.png');  
subplot(2,2,[1,3]);    imshow(I)
J=imrotate(I,35,'bilinear');
subplot(2,2,[2,4]);    imshow(J);

结果演示：

补充：图像存储格式使用函数imwrite.
$$imwrite(A,filename)$$
将图像数据 A 写入 filename 指定的文件，并从扩展名推断出文件格式。imwrite 在当前文件夹中创建新文件。输出图像的位深度取决于 A 的数据类型和文件格式。对于大多数格式来说：

• 如果 A 属于数据类型 uint8，则 imwrite 输出 8 位值。
• 如果 A 属于数据类型 uint16 且输出文件格式支持 16 位数据（JPEG、PNG 和 TIFF），则 imwrite 将输出 16 位的值。如果输出文件格式不支持 16 位数据，则 imwrite 返回错误。
• 如果 A 是灰度图像或者属于数据类型 double 或 single 的 RGB 彩色图像，则 imwrite 假设动态范围是 [0,1]，并在将其作为 8 位值写入文件之前自动按 255 缩放数据。如果 A 中的数据是 single，则在将其写入 GIF 或 TIFF 文件之前将 A 转换为 double。
• 如果 A 属于 logical 数据类型，则 imwrite 会假定数据为二值图像并将数据写入位深度为 1 的文件（如果格式允许）。BMP、PNG 或 TIFF 格式以输入数组形式接受二值图像。

如果 A 包含索引图像数据，则应另外指定 map 输入参数。

例：

imwrite(I,'pout2.png');

6.5使用MATLAB分析图像：目标计数

6.5.1图像预处理

$$Z=imsubtract(X,Y)$$

将数组Y中的每个元素从数组X中的相应元素中细分，并返回输出数组Z中相应元素中的差异。

如果X是整数数组，则输出中超出整数类型范围的元素将被截断，分数值将被舍入。

处理演示：

clc;clear;close all
I=imread('rice.png');
subplot(2,3,[1,4]);  imshow(I);
BG=imopen(I,strel('disk',15));
subplot(2,3,[2,5]);  imshow(BG);
I2=imsubtract(I,BG);
subplot(2,3,[3,6]);   imshow(I2)

演示结果：

6.5.2图像进行二值化处理

对比二值化前与二值化后之间的变化：

clc;clear;close all
I=imread('rice.png');   level=graythresh(I);
bw=im2bw(I,level);      subplot(2,2,[1,3]);
imshow(bw);            title('直接进行二值化')
BG=imopen(I,strel('disk',15));
I2=imsubtract(I,BG);   level=graythresh(I2);
bw2=im2bw(I2,level);
subplot(2,2,[2,4]);     imshow(bw2);
title('去除背景后进行二值化');

演示结果：

6.5.3目标计数:标记连通区域

识别米粒个数的关键在于识别连通区域.

MATLAB自带的bwlabel()函数计算连通区域,该函数使用了连通区域标记算法,将每个连通区域内的像素点赋值为同一个值.[简称区域生长]

$$L=bwlabel(BW)$$
返回LABEL MATRIX L，其中包含在BW中找到的8连接对象的标签。

使用GPU在二维二进制图像中标记连接的组件。

演示区域生长：

clc;clear;close all
I=imread('rice.png');
BG=imopen(I,strel('disk',15));
I2=imsubtract(I,BG);   level=graythresh(I2);
BW=im2bw(I2,level);
[labeled,numObjects]=bwlabel(BW,8);
numObjects   %一共有多少个

演示结果：

6.5.4以图片形式呈现

$$RGB=label2rgb(L)$$

将标签矩阵L（例如由标签矩阵、BWLABEL、BWLABELN或流域返回的矩阵L）转换为RGB颜色图像，以查看标签区域。LABEL2RGB函数根据标签矩阵中对象的数量确定要分配给每个对象的颜色。LABEL2RGB函数从颜色映射的整个范围中拾取颜色。

演示：

clc;clear;close all
I=imread('rice.png');
BG=imopen(I,strel('disk',15));
I2=imsubtract(I,BG);   level=graythresh(I2);
BW=im2bw(I2,level);
[labeled,numobjects]=bwlabel(BW,8);
RGB_label=label2rgb(labeled);  imshow(RGB_label);

演示结果：

[注：偏蓝色表示数值比较小，偏红色表示数值比较大]

6.5.5分析检测结果

$$stats=regionprops(BW,properties)$$

返回二进制图像BW中每个8连接组件（对象）的属性指定的属性集的度量值。stats是一个结构数组，其中包含图像中每个对象的结构

演示分析检测：

clc;clear;close all
I=imread('rice.png');
BG=imopen(I,strel('disk',15));
I2=imsubtract(I,BG);   level=graythresh(I2);
BW=im2bw(I2,level);
[labeled,numobjects]=bwlabel(BW,8);
graindata=regionprops(labeled,'basic');
graindata(51)

结果演示：

[注：演示第51颗米的面积，位置，边界框]

补充：选择连通区域
$$BW2=bwselect(BW,c,r,n)$$
返回包含重叠像素（r，c）的对象的二进制图像，其中n指定连接。对象是像素上的一组连接对象，即值为1的像素。默认情况下，b为四个连接的对象选择“查找”。

演示：

clc;clear;close all
I=imread('rice.png');
BG=imopen(I,strel('disk',15));
I2=imsubtract(I,BG);  BW=im2bw(I2,graythresh(I2));
objI=bwselect(BW);   imshow(objI);

演示结果：

7.数值微积分

7.1线性（多项式）数值运算

7.1.1使用向量表示多项式

在MATLAB中,多项式可以用向量表示,向量中的元素为多项式的系数(降幂排序):第一位为多项式最高次项系数,最后一位为常数项.

例如多项式:
$$f(x)=x^3-2x^2-5$$

可以用向量p = [1 0 -2 -5]表示.

7.1.2多项式求值

$$y=polyval(p,x)$$

计算多项式 p 在 x 的每个点处的值。参数 p 是长度为 n+1 的向量，其元素是 n 次多项式的系数（降幂排序）：

​ $$f(x)=f_1{x}^n+f_2{x}^{n-1}+...+f_{n}x+f_n+1$$

虽然可以为不同目的使用 polyint、polyder 和 polyfit 等函数计算 p 中的多项式系数，但您也可以为系数指定任何向量。

要以矩阵方式计算多项式，请改用 polyvalm。

演示：

clc;clear;close all
a=[9,-5,3,7];    x=-2:0.01:5;
f=polyval(a,x);
plot(x,f,'LineWidth',2);
xlabel('x');   ylabel('f(x)');
set(gca,'FontSize',14)

演示结果：

7.1.3多项式的微分

$$k=polyder(p)$$

返回 p 中的系数表示的多项式的导数，$$k(x)=\frac{d}{dx}p(x)$$.

演示把$$f(x)=5x^4-2x^2+1$$函数进行微分：

clc;clear;close all
p=[5 0 -2 0 1];
polyder(p)

演示结果：

得到的导数$$f^{\prime }(x)=20x^3-4x^2$$.

补充：用MATLAB对$$f(x)=5x^4-2x^2+1$$求$$f^{\prime }(7)$$?

clc;clear;close all
p=[5 0 -2 0 1];
polyval(polyder(p),7)

结果：

则：$$f^{\prime }(7)=6832$$.

7.1.4多项式积分

$$q=polyint(p,k)$$

使用积分常量 k 返回 p 中系数所表示的多项式积分。

演示对$$f(x)=5x^4-2x^2+1$$进行积分并求$$f(10)$$：

clc;clear;close all
p=[5 0 -2 0 1];
polyint(p)
polyval(polyint(p),7)

演示结果：

则积分$$f(x)=x^5-\frac{2}{3}{x}^3+x$$,$$f(10)=16585$$.

7.2非线性表达式的数值运算

7.2.1数值微分

7.2.1.1求差分

$$Y=diff(X)$$

计算沿大小不等于 1 的第一个数组维度的 X 相邻元素之间的差分：

• 如果 X 是长度为 m 的向量，则 Y = diff(X) 返回长度为 m-1 的向量。Y 的元素是 X 相邻元素之间的差分。

  $$Y=[x(2)-x(1)]\ast [x(3)-x(2)]\ast ...\ast [x(m)-x(m-1)]$$

• 如果 X 是不为空的非向量 p×m 矩阵，则 Y = diff(X) 返回大小为 (p-1)×m 的矩阵，其元素是 X的行之间的差分。

  $$Y=[x(2,:)-x(1,:);x(3,:)-x(2,:);...;x(p,:)-x(p-1,:)]$$

• 如果 X 是 0×0 的空矩阵，则 Y = diff(X) 返回 0×0 的空矩阵。

演示使用diff(X, n)计算向量X的n阶差分1,2,3,

clc;clear;close all
x=[1 2 5 2 1];
diff(x)   %n为1的时候
diff(x,2)   %n为2的时候
diff(x,3)   %n为3的时候

演示结果：

7.2.1.2求导数

使用导数的定义

​ $$f^{\prime }(x_0)=\underset{h\to 0}{\lim }\frac{f(x_0+h)-f(x_0)}{h}$$

可以计算函数在某点的近似导数.

演示给定$$f(x)=sin(x)$$，用$$h=0.1$$在$$x_0=\frac{\pi }{2}$$找到$$x_0$$:

clc;clear;close all
x0=pi/2;   h=0.1;
x=[x0 x0+h];
y=[sin(x0)   sin(x0+h)];
m=diff(y)./diff(x)

演示结果：

则数值为-0.05.

补充：h的取值越小,得到的近似导数越精确

下面程序计算$$f(x)=x^3$$的一阶和二阶导数的值并以图像的形式呈现出来.

clc;clear;close all
x=-2:0.005:2;  y=x.^3;
m=diff(y)./diff(x);
m2=diff(m)./diff(x(1:end-1));

plot(x,y,x(1:end-1),m,x(1:end-2),m2);
xlabel('x','FontSize',18);
ylabel('x','FontSize',18);
legend('f(x)=x^3','f','''''(x)',4);
set(gca,'FontSize',18)

演示结果：

7.2.2数值积分

7.2.2.1数值积分原理

有三种常见算法用于计算数值积分: 矩形法,梯形法,抛物线法,它们分别把微分区间的图形视为矩形,梯形,抛物线以计算面积.

算法	图示	表达式
矩形法(Midpoint Rule)		$$f_a^{b}f(x)dx=2h\sum _{i=0}^{(n/2)-1}f(x_{2i+1})$$
梯形法(Trapezoid Rule)		$$f_a^{b}f(x)dx=\frac{h}{2}[f(a)+2\sum _{i=1}^{n-1}f(x_i)+f(b)]$$
抛物线法(Simpson’s Rule)		$$f_a^{b}f(x)dx=\frac{h}{3}[f(a)+2\sum _{i=1}^{n/2-1}f(x_{2i})+4\sum _{i=1}^{n/2}f(x_{2}i-1)+f(b)]$$

7.2.2.2矩形法计算

使用和的中点规则

$${\int }_{x_0}^{x_3}f(x)dx\approx h{f}_0+h{f}_1+h{f}_2=h\sum _{i=1}^{n-1}{f}_i$$

[关键函数：midpoint]

例：求$$A={\int }_0^{2}4x^{3}dx$$.

思路：在0—2之间设定每0.05为一个矩形面积，找出中心点求出所有单个面积，再进行求和。

clc;clear
h=0.05;    x=0:h:2;
midpoint=(x(1:end-1)+x(2:end))./2;
y=4*midpoint.^3;
s=sum(h*y)

演示结果：

则精确结果为15.995，趋近于16

7.2.2.3梯形计算法

$${\int }_{x_0}^{x_3}f(x)dx\approx h\frac{f_0+f_1}{2}+h\frac{f_1+f_2}{2}+h\frac{f_2+f_3}{2}=h\sum _{i=0}^{n-1}fi$$

[关键函数：trapz]

例：求$$A={\int }_0^{2}4x^{3}dx$$.

思路：分出有限个梯形，使用上地加下地乘以高除以2，求出各个梯形的面积，然后把每一个梯形面积相加。

clc;clear
h=0.05;   x=0:h:2;
y=4*x.^3;
s=h*trapz(y)

演示结果：

则精确结果为16.0100，趋近于16

7.2.2.4抛物线计算法

$${\int }_{x_0}^{x_2}f(x)dx\approx \frac{h}{3}(f_0+4f_1+f_2)$$.

[关键函数：Rule]

例：求$$A={\int }_0^{2}4x^{3}dx$$.

思路：把整个面积平均分成两部分，起始位置的高度加四倍分裂位置的高度加末端位置的高度，三者相加乘以第一部分的长度，再除以3，取平均面积。

clc;clear
h=0.05;   x=0:h:2;
y=4*x.^3;
s=h/3*(y(1)+2*sum(y(3:2:end-2))+4*sum(y(2:2:end))+y(end))

演示结果：

则精确结果为16

总结：抛物线计算法的精准度比矩形计算法和梯形计算法更加准确。

7.3查看函数句柄（@）

句柄是指向函数的指针，可用于将函数传递给其他函数。

例如，以下函数的输入是另一个函数

function [y]=x_plot(input,x)
% xy plot recciyos the handle of a function and plots that
% function of x
y=input(x);   plot(x,y,'r--');
xlabel('x');   ylabel('function(x)');
end

在行窗命令口输入：

x_plot(@sin,0:0.01:2*pi)

x_plot(@cos,0:0.01:2*pi)

结果演示：

补充：
$$q=integral(fun,xmin,xmax)$$
使用全局自适应积分和默认误差容限在 xmin 至 xmax 间以数值形式为函数 fun 求积分。

1.练习 $${\int }_0^2\frac{1}{x^3-2x-5}dx$$.

y = @(x) 1./(x.^3-2*x-5);
integral(y,0,2)

演示结果： -0.4605

2.练习 $${\int }_0^{\pi }{\int }_{\pi }^{2\pi }(y\ast sin(x)+x\ast cos(y))dxdy$$.

f = @(x,y) y.*sin(x)+x.*cos(y);
integral2(f,pi,2*pi,0,pi)

演示结果： -9.8696

3.练习 $${\int }_{-1}^1{\int }_0^1{\int }_0^{\pi }(y\ast sin(x)+z\ast cos(y))dxdydz$$.

f = @(x,y,z) y.*sin(x)+z.*cos(y);
integral3(f,0,pi,0,1,-1,1)

演示结果： 2.0000

8.方程式求根

8.1创建符号数字

$$symsvar1...varN$$

创建符号变量var1。。。。用空格分隔不同的变量。syms从变量中清除所有假设。

演示：

syms x
x+x+x
(x+x+x)/4

演示结果： $$3x$$和$$\frac{3x}{4}$$.

8.1.1符号寻根

$$S=solve(eqn,var)$$

求解变量var的方程eqn。如果未指定var，symvar函数将确定要求解的变量。例如，solve（x+1==2，x）为x解方程x+1=2。

演示求$$y=x\ast sin(x)-x=0$$求解：

clc;clear
syms x 
y=x*sin(x)-x;
solve(y,x)

演示结果：$$x=0$$或$$\frac{\pi }{2}$$.

8.1.2解多重方程

例1：$$\{\begin{array}{l}x-2y=5\\ x+y=6\end{array}$$.

clc;clear
syms x y
eq1=x-2*y-5;
eq2=x+y-6;
A=solve(eq1,eq2,x,y);
A.x   
A.y

演示结果：$$x=\frac{17}{3}$$ 和$$y=\frac{1}{3}$$.

8.1.3解方程零数组

$$x=fsolve(fun,x0)$$

从x0开始，尝试求解方程fun（x）=0，一个零数组。

演示$$f(x)=1.2x+0.3+x\ast sin(x)$$：

clc;clear
f2 = @(x) (1.2*x+0.3+x*sin(x));
fsolve(f2,0)

演示结果： -0.3500

8.1.4多项式求根

$$r=roots(p)$$

以列向量的形式返回 p 表示的多项式的根。输入 p 是一个包含 n+1 多项式系数的向量，以 xn 系数开头。0 系数表示方程中不存在的中间幂。例如：p = [3 2 -2] 代表多项式 3x2+2x−2。

roots 函数对 $$p_{\ast }{x}^n+...+p_{n}x+p_{n+1}=0$$格式的多项式方程求解。包含带有非负指数的单一变量的多项式方程。

演示求$$f(x)=x^5-3.5x^4+2.75x^3+2.125x^2-3.875x+1.25$$:

clc;clear
roots([1 -3.5 2.75 2.125 -3.875 1.25])

演示结果：

8.1.5符号表达式的化简与代入

使用simplify()函数可以化简符号表达式.

syms x a b c

simplify(sin(x)^2 + cos(x)^2);		% 得到 1
simplify(exp(c*log(sqrt(a+b))));	% 得到 (a + b)^(c/2)
1234

表达式化简的标准是不确定的,下面三个函数分别按照不同标准化简表达式:

• expand()函数可以展开表达式

  syms x
  
  f = (x ^2- 1)*(x^4 + x^3 + x^2 + x + 1)*(x^4 - x^3 + x^2 - x + 1);
  expand(f);		% 得到 x^10 - 1
  1234
  

• factor()函数可以分解因式

  syms x
  
  g = x^3 + 6*x^2 + 11*x + 6;
  factor(g);		% 得到 (x + 3)*(x + 2)*(x + 1)
  1234
  

• horner()函数可以将多项式变为嵌套形式

  syms x
  
  h = x^5 + x^4 + x^3 + x^2 + x;
  horner(h);		% 得到 x*(x*(x*(x*(x + 1) + 1) + 1) + 1)
  1234
  

8.1.6符号表达式的代入

使用sub(expr, old, new)函数可以将符号表达式expr中的old替换为new.

syms x

f = 2*x^2 - 3*x + 1;
subs(f, 1/3)	% 得到 2/9
1234
syms x y

f = x^2*y + 5*x*sqrt(y);
subs(f, x, 3);	% 得到 9*y + 15*y^(1/2)

8.1.7符号微积分运算

8.1.7.1极限

使用limit(expr, var, a)函数可以求取符号表达式expr在变量var趋近于a时的极限,添加参数'left'或'right'可以指定左极限或右极限.

syms x;
expr = 1/x;

limit(expr,x,0);			% 得到NaN
limit(expr,x,0,'left');		% 得到-Inf
limit(expr,x,0,'right');	% 得到Inf
123456

8.1.7.2微分

使用diff(expr, var, n)函数可以求取符号表达式expr对变量var的n阶微分.

syms a b c x;
expr = a*x^2 + b*x + c;

diff(expr, a);		% 得到 x^2
diff(expr, b);		% 得到 x
diff(expr, x);		% 得到 b + 2*a*x
diff(expr, x, 2);	% 得到 2*a
1234567

8.1.7.3积分

使用int(expr, var)函数可以求取符号表达式expr对变量var的不定积分.使用int(expr, var, [a, b])函数可以指定上下限求定积分,a和b可以是符号表达式.

syms x a b
expr = -2*x/(1+x^2)^2;

int(expr, x);			% 得到 1/(x^2 + 1)
int(expr, x, [1, 2]);	% 得到 -0.3
int(expr, x, [1, Inf]);	% 得到 -0.5
int(expr, x, [a, b]);	% 得到 1/(b^2 + 1) - 1/(a^2 + 1)
1234567

对于一些函数,MATLAB不能求出其积分,这时MATLAB会返回一个未解析(unsolved)的积分形式.

syms x
int(sin(sinh(x)));	% 一个无解的积分,MATLAB返回 int(sin(sinh(x)), x) 
12

8.1.7.4级数求和

使用symsum(expr, k, [a b])计算级数expr的索引k从a到b的加和.

syms k x

symsum(k^2, k)			% 得到 k^3/3 - k^2/2 + k/6
symsum(k^2, k, [0 10])	% 得到 385
symsum(x^k/factorial(k),k,1,Inf)	% 得到 exp(x) - 1
12345

8.1.7.5泰勒展开

使用taylor(expr,var,a)计算表达式expr在var=a处的泰勒级数.

syms x

taylor(exp(x))	% 得到 x^5/120 + x^4/24 + x^3/6 + x^2/2 + x + 1
taylor(sin(x))	% 得到 x^5/120 - x^3/6 + x
taylor(cos(x))	% 得到 x^4/24 - x^2/2 + 1

8.2二等分法

通过迭代方法改变范围，确定位置用逻辑用于来说：

8.3牛顿迭代法

牛顿迭代法的一个迭代关系式：

$$x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)}$$

如果是连续的，并且待求的零点是孤立的，那么在零点周围存在一个区域，只要初始值位于这个邻近区域内，那么牛顿法必定收敛。 并且，如果不为0, 那么牛顿法将具有平方收敛的性能. 粗略的说，这意味着每迭代一次，牛顿法结果的有效数字将增加一倍。

迭代法，是一种不断用变量的旧值递推新值的过程，跟迭代法相对应的是直接法（或者称为一次解法），即一次性解决问题。迭代算法是用计算机解决问题的一种基本方法。它利用计算机运算速度快、适合做重复性操作的特点，让计算机对一组指令（或一定步骤）重复执行，在每次执行这组指令（或这些步骤）时，都从变量的原值推出它的一个新值。

8.4绘制图像

可以对符号表达式绘制图像,常用的绘图函数如下:

函数	作用
fplot()	绘制符号表达式的二维线图像
fplot3()	绘制符号表达式的三维线图像
ezpolar()	绘制符号表达式的极坐标线图像
fmesh()	绘制网状面图像
fsurf()	绘制带颜色的面图像
fcontour()	绘制轮廓图像
fimplicit()	绘制隐含函数关系的图像

下面例子展示二维和三维线图像的绘制

clc;clear;close all
subplot(1, 2, 1)
syms x
f = x^3 - 6*x^2 + 11*x - 6;
fplot(f, x)

subplot(1, 2, 2)
syms t
fplot3(t^2*sin(10*t), t^2*cos(10*t), t)

下面例子演示三维面的绘制:

clc;clear;close all
syms x y
fsurf(x^2 + y^2)

下面例子演示隐含函数关系图像的绘制:

clc;clear;close all
syms x y
eqn = (x^2 + y^2)^4 == (x^2 - y^2)^2;
fimplicit(eqn, [-1 1])

9.线性规划与线性系统

9.1高斯消去法

假设给定：

$$\{\begin{array}{l}x+2y+z=2\\ 2x+6y+z=7\\ x+y+4z=3\end{array}$$. $$\Rightarrow$$ $$\left[\begin{array}{cccc}1& 2& 1& 2\\ 2& 6& 1& 7\\ 1& 1& 4& 3\end{array}\right]$$.

$$\Rightarrow$$ $$\left[\begin{array}{cccc}1& 2& 1& 2\\ 0& 2& -1& 3\\ 0& -1& 3& 1\end{array}\right]$$.$$\Rightarrow$$ $$\left[\begin{array}{cccc}1& 2& 1& 2\\ 0& 2& -1& 3\\ 0& 0& 5/2& 5/2\end{array}\right]$$

高斯消去法，也称为逐次淘汰。

9.1.1MATLAB运算高斯消去法

$$R=rref(A)$$

使用 Gauss-Jordan 消元法和部分主元消元法返回简化行阶梯形的 A。

演示：$$\{\begin{array}{l}x+2y+z=2\\ 2x+6y+z=7\\ x+y+4z=3\end{array}$$. $$\Rightarrow$$ $$\left[\begin{array}{cccc}1& 2& 1& 2\\ 2& 6& 1& 7\\ 1& 1& 4& 3\end{array}\right]$$.

clc;clear;close all
A = [1 2 1;2 6 1;1 1 4];
B = [2 7 3]';
R = rref([A B])

演示结果：

 1     0     0    -3
 0     1     0     2
 0     0     1     1

则结果化为标准形式为：$$\left[\begin{array}{cccc}1& 0& 0& -3\\ 0& 1& 0& 2\\ 0& 0& 1& 1\end{array}\right]$$.

9.2LU因子分解

假设我们想求解：$$Ax=b$$，其中$$A\in {\rho }^{m\times m}$$.

将A分解为2个三角形矩阵：$$A=L^{-1}U$$.

问题变为：$$AX=B\Rightarrow L^{-1}Ux=b$$.

策略

1.求解$$L^{-1}y=b$$得到$$y$$.

2.然后解出$$Ux=y$$.

LU现实化实例:$$A=\left[\begin{array}{ccc}1& 1& 1\\ 2& 3& 5\\ 4& 6& 8\end{array}\right]$$.

$$L_{1}A=\left[\begin{array}{ccc}1& 0& 0\\ -2& 1& 0\\ -4& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 1& 1\\ 2& 3& 5\\ 4& 6& 8\end{array}\right]=\left[\begin{array}{ccc}1& 1& 1\\ 0& 1& 3\\ 0& 2& 4\end{array}\right]$$.

$$L_2(L_{1}A)=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& -2& 1\end{array}\right]\left[\begin{array}{ccc}1& 1& 1\\ 0& 1& 3\\ 0& 2& 4\end{array}\right]=\left[\begin{array}{ccc}1& 1& 1\\ 0& 1& 3\\ 0& 0& -2\end{array}\right]=U$$

9.2.1MATLAB运输LU因式分解

$$[L,U]=lu(A)$$

将满矩阵或稀疏矩阵 A 分解为一个上三角矩阵 U 和一个经过置换的下三角矩阵 L，使得 A = L*U。

演示：$$A=\left[\begin{array}{ccc}1& 1& 1\\ 2& 3& 5\\ 4& 6& 8\end{array}\right]$$.

clc;clear;close all
A=[1 1 1;2 3 5;4 6 8];
[L,U,P]=lu(A);
inv(L)
U

演示结果：

4.0000    6.0000    8.0000
     0   -0.5000   -1.0000
     0         0    1.0000

则结果：$$\left[\begin{array}{ccc}4& 6& 8\\ 0& -0.5& -1\\ 0& 0& 1\end{array}\right]x=y$$.

9.3矩阵左除：\or mldivide()

使用因子分解方法求解线性方程组$$Ax=b$$：

$$\{\begin{array}{l}x+2y+z=2\\ 2x+6y+z=7\\ x+y+4z=3\end{array}$$.

演示：

clc;clear;close all
A=[1 2 1;2 6 1;1 1 4];
b=[2;7;3];
x=A\b

演示结果：

-3.0000
 2.0000
 1.0000

则结果：$$x=[-3,2,1{]}^{\prime }$$.

9.4矩阵分解函数

代码	解析
qr	正交三角分解
ldl	Hermitian不定矩阵的块LDL分解
ilu	稀疏不完全因子分解
chol	Cholesky因子分解
gsvd	广义奇异balue分解
svd	单值分解

克莱默（逆）法：通过矩阵的秩来逆矩阵。

给定方程：$$\{\begin{array}{l}x+2y+z=2\\ 2x+6y+z=7\\ x+y+4z=3\end{array}⟹\left[\begin{array}{ccc}1& 2& 1\\ 2& 6& 1\\ 1& 1& 4\end{array}\right]x=\left[\begin{array}{c}2\\ 7\\ 3\end{array}\right]$$

$$x=A^{-1}b$$

演示：

clc;clear;close all
A=[1 2 1;2 6 1;1 1 4];
b=[2;7;3];
x=inv(A)*b

演示结果：

   -3.0000
    2.0000
    1.0000

则结果为：$$x=\left[\begin{array}{c}-3\\ 2\\ 1\end{array}\right]$$.

补充：逆矩阵不是无时无刻都存在，通常有3种情况，1.存在且唯一2.存在且无数解3.不存在

9.5线性系统

假设你得到了线性方程组:

$$\{\begin{array}{l}2\cdot 2-12\cdot 4=x\\ 1\cdot 2-5\cdot 4=y\end{array}$$.

矩阵表示法：

$$\left[\begin{array}{cc}2& -12\\ 1& -5\end{array}\right]\left[\begin{array}{c}2\\ 4\end{array}\right]=\left[\begin{array}{c}x\\ y\end{array}\right]$$.

特征值与特征向量

对于系统$$A\in {\rho }^{m\times m}$$，矩阵乘法$$y=Ab$$是复杂的

想找到向量$$v_i\in {\rho }^m$$.

$$A{v}_i={\lambda }_i{v}_i$$， 哪里$${\lambda }_i\in {\rho }^m$$

然后分解$$b=\sum {\alpha }_i{v}_i$$,$${\alpha }_i\in \rho$$.

复制变成:$$Ab=\sum {\alpha }_{i}A{v}_i=\sum {\alpha }_i{\lambda }_i{v}_i$$.

求解特征值和特征向量:

例:$$Ab=\left[\begin{array}{cc}2& -12\\ 1& -5\end{array}\right]\left[\begin{array}{c}2\\ 4\end{array}\right]$$.

$${\lambda }_1{v}_1=-1\left[\begin{array}{c}0.97\\ 0.24\end{array}\right]$$, $${\lambda }_2{v}_2=-2\left[\begin{array}{c}0.95\\ 0.32\end{array}\right]$$,

$$b={\alpha }_1{v}_1+{\alpha }_2{v}_2=-41.2v_1+44v_2$$.

$$Ab=A({\alpha }_1{v}_1+{\alpha }_2{v}_2)={\alpha }_{1}A{v}_1+{\alpha }_{2}A{v}_2={\alpha }_1{\lambda }_1{v}_1+\alpha {\lambda }_2{v}_2=(-41.2)(-1)\left[\begin{array}{c}0.97\\ 0.24\end{array}\right]+(44.3)(-2){\lambda }_2{v}_2=-2\left[\begin{array}{c}0.95\\ 0.32\end{array}\right]$$

9.5.1MATLAB进行线性系统计算

$$[V,D]=eig(A,B)$$

返回广义特征值的对角矩阵 D 和满矩阵 V，其列是对应的右特征向量，使得 A*V = B*V*D。

演示$$A=\left[\begin{array}{cc}2& -12\\ 1& -5\end{array}\right]$$：

clc;clear
A=[2 -12;1 -5];
[v,d]=eig(A)

演示结果：

v =

    0.9701    0.9487
    0.2425    0.3162

d =

   -1.0000         0
         0   -2.0000

则结果特征值：$$\left[\begin{array}{cc}0.9701& 0.9487\\ 0.2425& 0.3162\end{array}\right]$$,特征向量：$$\left[\begin{array}{cc}-1& 0\\ 0& -2\end{array}\right]$$.

9.6矩阵指数表达式

$$Y=expm(X)$$

计算 X 的矩阵指数。虽然不按此种方式计算，但是如果 X 包含一组完整的特征向量 V 和对应特征值 D，则 [V,D] = eig(X) 且

expm(X) = V*diag(exp(diag(D)))/V

一个典型的线性时不变系统通常是制定为:

$$y=\frac{dx(t)}{dt}=x=Ax$$.

例$$A=\left[\begin{array}{ccc}0& -6& -1\\ 6& 2& -16\\ -5& 20& -10\end{array}\right]$$以3D图像的形式呈现：

clc;clear
A=[0 -6 -1;6 2 -16;-5 20 -10];
x0=[1 1 1]' ； X=[]‘;
for t=0:.01:1
    X=[X expm(t*A)*x0];
end
plot3(X(1,:),X(2,:),X(3,:),'-o');
xlabel('x_1');  ylabel('x_2');
zlabel('x_3');   grid on;
axis tight square

演示结果：

10.统计

科学统计包括数据的收集、分析、解释、呈现和组织。主要的统计方法：1.统计2.描述性统计3.推断统计

10.1集中趋势

MATLAB命令	解析
mean	数值平均值
median	数组中值
mode	数组中最频繁的值
prctile	数据集的百分位数
max	数组中最大的元素
min	数组中最小的元素
std	标准差
var	方差

练习使用上述命令：

clc;clear
load stockreturns
x4 = stocks(:,4);
mean(x4);
median(x4);
mode(x4)
prctile(x4,[1 3])
max(x4)
min(x4)
std(x4)
var(x4)

求得结果：

ans =
  -5.8728e-04

ans =
    0.0617

ans =
   -5.8764

ans =
   -5.6491   -4.6689

ans =
    6.4692

ans =
   -5.8764

ans =
    2.3590

ans =
    5.5649

数字总是更有力的,例：

x=1:14;
freqy = [1 0 1 0 4 0 1 0 3 1 0 0 1 1];
subplot(2,3,[1 4]);   bar(x,freqy); xlim([0 15]);
subplot(2,3,[2 5]);   area(x,freqy); xlim([0 15]);
subplot(2,3,[3 6]);   stem(x,freqy); xlim([0 15]);

箱线图可以突出显示数据的四分位点.

clc;clear
marks = [80 81 81 84 88 92 92 94 96 97];
boxplot(marks)
prctile(marks, [25 50 75])	% 得到 [81 90 94]

演示结果：

10.2变化情况

10.2.偏度

$$y=skewness(X)$$

偏度反映数据的对称程度

• 数据左偏时,其偏度小于0.
• 数据完全对称时,其偏度等于0.
• 数据右偏时,其偏度大于0.

演示：

clc;clear;close all
X=randn([10 3])*3;
X(X(:,1)<0,1)=0; X(X(:,3)>0,3)=0;
boxplot(X,{'Right-skewed','Symmetric','Left-skewed'});
y=skewness(X)

演示结果：

y =

0.3901   -0.8070   -0.8618

【注：每一次演算的结果不一所画出来的图也不一样】

10.3峰度

$$k=kurtosis(X)$$

返回X的样本峰度。

如果X是向量，则峰度（X）返回一个标量值，即X中元素的峰度。

如果X是矩阵，则峰度（X）返回一个行向量，该行向量包含X中每列的样本峰度。

如果X是一个多维数组，则峰度（X）沿X的第一个非单量子维数进行运算。

峰度(Kurtosis)表征概率密度分布曲线在平均值处峰值的高低.直观来看,峰度反映了峰部的尖度.

10.3 统计推断

推论统计的核心即为假设检验.下列函数用于进行假设检验.

函数	作用
ttest()	进行T检验
ztest()	进行Z检验
ranksum()	进行秩和检验
signrank()	进行符号秩检验

确定假设检验的可能性，比如说0.95

找到H的95%置信区间

检查你的分数是否在区间内

演示：两种股票收益率的平均值相同吗？

clc;clear;close all
load stockreturns 
x1=stocks(:,3);   x2=stocks(:,10);
boxplot([x1,x2],{'3','10'});
[h,p]=ttest2(x1,x2)

演示结果：

h =
     1

p =
    0.0423

【注：p和h我们一般的使用的范围都是在0.5—1，这里的p小于0.5，我们不使用，表示在默认显著性水平(5%)以上拒绝】

补充：

load examgrades
x = grades(:,1);
y = grades(:,2);
[h,p] = ttest(x,y);

结果：

h =
     0

p =
    0.9805

执行上述程序,得到[h p] = [0 0.9805],表示在默认显著性水平(5%)下我们没有理由拒绝x与y同分布.

11.回归与内插

11.1简单线性回归

简单线性回归为一元线性回归模型，是指模型中只含有一个自变量和因变量，该模型的数学公式可以表示为y=ax+b+ε，a为模型的斜率，b为模型的截距，ε为误差。

收集了一组数据点（x，y）

假设x和x是线性相关的

定义误差平方和

$$SSE=\sum _i{\epsilon }_i^2=\sum _i(y_i-\widehat{y_i}{)}^2$$

假设回归模型:$$\widehat{y_i}={\beta }_0+{\beta }_{1}x$$,

$$SSE=\sum _i(y_i-{\beta }_0-{\beta }_{1}x{)}^2$$.

当SSE相对于每个参数的梯度等于零时，SSE最小化:

$$\{\begin{array}{l}\frac{\partial {\sum }_i{ϵ}^2}{\partial {\beta }_0}=-2\sum _i(y_i-{\beta }_0-{\beta }_1{x}_1)=0\\ \frac{\partial {\sum }_i{ϵ}^2}{\partial {\beta }_1}=-2\sum _i(y_i-{\beta }_0-{\beta }_1{x}_1)x_1=0\end{array}$$.联立求解

假设存在N个数据点：

$$\{\begin{array}{l}\sum _{i=1}^N{y}_i={\beta }_0\cdot N+{\beta }_1\sum _{i=1}^N{x}_i\\ \sum _{i=1}^N{y}_i{x}_i={\beta }_0\sum _{i=1}^N{x}_i+{\beta }_1\sum _{i=1}^N{x}_i^2\end{array}⟹\left[\begin{array}{c}\sum y_i\\ \sum y_i{x}_i\end{array}\right]=\left[\begin{array}{cc}N& \sum x_i\\ \sum x_i& \sum x_i^2\end{array}\right]\left[\begin{array}{c}{\beta }_0\\ {\beta }_1\end{array}\right]$$

11.2一元多项式拟合

$$p=polyfit(x,y,n)$$

返回次数为 n 的多项式 p(x) 的系数，该阶数是 y 中数据的最佳拟合（在最小二乘方式中）。p 中的系数按降幂排列，p 的长度为 n+1

$$f(x)=p_1{x}^n+p_2{x}^{n-1}+...+p_{n}x+p_{n+1}$$.

演示使用polyfit(x, y, n)函数对数据x和y进行n次多项式拟合：

clc;clear;close all
x = [-1.2 -0.5 0.3 0.9 1.8 2.6 3.0 3.5];
y = [-15.6 -8.5 2.2 4.5 6.6 8.2 8.9 10.0];
fit=polyfit(x,y,1)
xfit=[x(1):0.1:x(end)];   yfit=fit(1)*xfit+fit(2);
plot(x,y,'ro',xfit,yfit);
set(gca,'FontSize',14);
legend('data pointa','best-fit');

演示结果：

fit =

    4.9897   -4.4491

补充：要精确一个线性是否准确需要多次线性来测量，引出万能线性代码。

clc;clear;close all
x = [-1.2 -0.5 0.3 0.9 1.8 2.6 3.0 3.5];
y = [-15.6 -8.5 2.2 4.5 6.6 8.2 8.9 10.0];
for i=1:3
    p = polyfit(x,y,i);	% 分别进行一次,二次,三次拟合
    xfit = x(1):0.1:x(end); yfit = polyval(p,xfit);
    subplot(1,3,i); plot(x,y,'ro',xfit,yfit); 
    legend('Data points','Fitted curve', 'Location', 'southeast');
end

演示结果3次拟合图：

显然第三次线性拟合效果最好，不同图像需要的线性程度不同。

11.3多元线性拟合

$$b=regress(y,X)$$

返回矩阵X中预测项向量y中响应的多元线性回归的系数估计向量b。若要计算具有常数项（截距）的模型的系数估计，请在矩阵X中包含一列系数估计。

演示：使用regress(y, X)函数对数据X和y进行多元线性回归.

clc;clear;close all
load carsmall;
y=MPG ; x1=Weight;  x2=Horsepower;     %导入数据集
X=[ones(length(x1),1) x1 x2];          %构建增广矩阵
b=regress(y,X);                        %进行线性回归
x1fit=min(x1):100:max(x1);             %开始绘图
x2fit=min(x2):100:max(x2);
[X1fit,X2fit]=meshgrid(x1fit,x2fit);
YFIT=b(1)+b(2)*X1fit+b(3)*X2fit;
scatter3(x1,x2,y,'filled');
hold on
mesh(X1fit,X2fit,YFIT);
hold off
xlabel('Weight');  ylabel('Horsepowed');   zlabel('MPG');
view(50,10);                         %调整显示角度

演示结果：

11.4非线性拟合

对于非线性拟合,需要使用曲线拟合工具箱.在命令窗口输入cftool()打开曲线拟合工具箱.

11.5一维插值

下列函数均与一维插值有关:

函数	作用
interp1(x,v)或interp1(x,v,xq)	线性插值
spline(x,v)或spline(x,v,xq)	三次样条插值
pchip(x,v)或pchip(x,v,xq)	三次Hermite插值
mkpp(breaks,coefs)	生成分段多项式
ppval(pp,xq)	计算分段多项式的插值结果

下面例子演示使用interp1(x, v, xq)进行线性插值和使用spline(x, v, xq)进行三次样条插值.各参数意义如下:

• x,v: 待插值样本点.
• xq: 查询点,函数返回在这些点处的插值结果.

clc;clear;close all
% 构造数据
x = linspace(0, 2*pi, 40); x_m = x; x_m([11:13, 28:30]) = NaN; 
y_m = sin(x_m);
plot(x_m, y_m, 'ro', 'MarkerFaceColor', 'r'); hold on; 	

% 对数据进行线性插值
m_i = ~isnan(x_m);
y_i = interp1(x_m(m_i), y_m(m_i), x);
plot(x,y_i, '-b'); hold on;

% 对数据进行三次样条插值
m_i = ~isnan(x_m);
y_i = spline(x_m(m_i), y_m(m_i), x);
plot(x,y_i, '-g'); 
legend('Original', 'Linear', 'Spline');

演示结果：

三次样条插值的原理在每两个样本点之间用两两相切的三次函数曲线来插值.若spline()方法不指定查询点,则返回一个结构体,其中封装了插值三次函数的系数.使用mkpp(breaks, coefs)可以创建一个类似的结构体,使用ppval(pp, xq)可以计算查询点的插值结果.

使用pchip(x, y, xq)函数可以进行三次Hermite插值,该算法同样以三次函数进行插值,但得到的曲线更平缓.

clc;clear;close all
x = -3:3; y = [-1 -1 -1 0 1 1 1]; xq1 = -3:.01:3;
p = pchip(x,y,xq1);
s = spline(x,y,xq1);
plot(x,y,'o',xq1,p,'-',xq1,s,'-.')
legend('Sample Points','pchip','spline','Location','SouthEast')

演示结果：

11.6二维插值

使用interp2()可以进行二维插值,向其method参数传入字符串可以指定插值算法.

方法	说明
‘linear’	(默认)在查询点插入的值基于各维中邻点网格点处数值的线性插值.
‘spline’	在查询点插入的值基于各维中邻点网格点处数值的三次插值.插值基于使用非结终止条件的三次样条.
‘nearest’	在查询点插入的值是距样本网格点最近的值.
‘cubic’	在查询点插入的值基于各维中邻点网格点处数值的三次插值.插值基于三次卷积.
‘makima’	修改后的Akima三次Hermite插值.在查询点插入的值基于次数最大为3的多项式的分段函数,使用各维中相邻网格点的值进行计算.为防过冲,已改进 Akima 公式.

clc;clear;close all
% 构建样本点
xx = -2:.5:2; yy = -2:.5:3; [x,y] = meshgrid(xx,yy); 
xx_i = -2:.1:2; yy_i = -2:.1:3; [x_i,y_i] = meshgrid(xx_i,yy_i);
z = x.*exp(-x.^2-y.^2);

% 线性插值
subplot(2, 2, [1 3]); 
z_i = interp2(xx,yy,z,x_i,y_i);
surf(x_i,y_i,z_i); hold on;
plot3(x,y,z+0.01,'ok','MarkerFaceColor','r'); hold on;

% 三次插值
subplot(2, 2, [2 4]); 
z_ic = interp2(xx,yy,z,x_i,y_i, 'spline');
surf(x_i,y_i,z_ic); hold on;
plot3(x,y,z+0.01,'ok','MarkerFaceColor','r'); hold on;