matlab语言与应用 09 概率论与数理统计

现代科学运算-matlab语言与应用

$\qquad \qquad$ 东北大学 http://www.icourses.cn/home/ (薛定宇)
《高等应用数学问题的MATLAB求解（第四版）》
代码可在matlab r2016b 运行。

09 概率论与数理统计问题的计算机求解

09.01 常用概率分布

概率与概率密度函数
什么是概率？
随机事件发生的可能性
如果随机变量$\xi$落入任意(a, b)区间的概率
$F(a \leq \xi \leq b) = \int_a^b p(x) dx$
则称 $p(x)$ 为概率密度函数, 满足
$p(x) \geq 0, \int_{-\infty}^{\infty} p(x) dx = 1$

概率分布函数
由概率密度可以定义概率分布函数
$F(x) = \int_{-\infty}^{x} p(t) dt$
物理意义：$概率分布函数F(x) - 随机变量 \xi 满足 \xi \leq x 发生的概率$
分布函数性质：
函数F(x)单调非减函数
$0 \leq F(x) \leq 1, F(-\infty) = 0, F(\infty) = 1$

常见分布的概率密度函数与分布函数
Possion分布
征态分布
F分布
T分布
$\chi^2$分布
Gamma分布
Rayleigh分布

后缀: pdf, cdf, inv, rnd, stat, fit
例如：batapdf(), gamcdf(), raylrnd()
pdf(‘beta’, …), cdf(‘gam’, …), random(‘rayl’, …)
inv->icdf, rnd->random, fit->fittest
$\begin{array}{l l l | l l l | l l l} \hline 关键词 & 分布名称 & 有关参数 &关键词 & 分布名称 & 有关参数 & 关键词 & 分布名称 & 有关参数 \\ \hline beta & Beta分布 & a, b & bino & 二项式分布 & n, p & chi2 & \chi^2分布 & k \\ ev & 极值分布 & \mu, \sigma & exp & 指数分布 & \lambda & f & F分布 & p, q \\ gam & Gamma分布 & a, \lambda & geo & 几何分布 & p & hyge & 超几何分布 & m, p, n \\ logn & 对数征态分布 & \mu, \sigma & mvn & 多变量征态分布 & \mu, \sigma & nbin & 负二项式分布 & v_1, v_2, \delta \\ ncf & 非零F分布 & k, \delta & nct & 非零T分布 & k, \delta & ncx2 & 非零\chi^2分布 & k, \delta \\ norm & 征态分布 & \mu, \sigma & poiss & Poisson分布 & \lambda & ray1 & Rayleigh分布 & b \\ t & T分布 & k & unif & 均匀分布 & a, b & wbl & Webull分布 & a, b \\ \hline \end{array}$

Poisson分布
Poisson分布的概率密度为:
$p(x) = \dfrac{\lambda^x}{x!} e^{-\lambda x}, x = 0, 1, 2, 3, \cdots$
其中, $\lambda$ 为正整数
Poisson分布的概率密度函数:
y = poisspdf(x, $\lambda$), y = pdf(‘poiss’, x, $\lambda$)
F = poisscdf(x, $\lambda$), F = cdf(‘poiss’, x, $\lambda$)
x = poissinv(F, $\lambda$), x = icdf(‘poiss’, F, $\lambda$)

例9-1 Poisson分布
Poisson分布的概率密度函数与分布函数曲线
Poisson分布的参数 $\lambda = 1, 2, 5, 10$

x = [0:15]';
y1 = [];
y2 = [];
lam1 = [1, 2, 5, 10];
for i = 1:length(lam1)
    y1 = [y1, poisspdf(x, lam1(i))];
    y2 = [y2, poisscdf(x, lam1(i))];
end;
stem(x, y1);
line(x, y1);
figure;
stem(x, y2);
line(x, y2);

也可以使用pdf()、cdf()、icdf()

征态分布
征态分布的概率密度函数为:
$p(x) = \dfrac{1}{\sqrt{2 \pi}\sigma} e^{-\dfrac{(x-\mu)^2}{2\sigma^2}}$
其中，$\mu 和 \sigma^2$ 分别为征态分布的均值和方差
征态分布的概率密度函数调用格式：
$y = normpdf(x, \mu, \sigma) \quad y = pdf('norm', x, \mu, \sigma)$
$F = normcdf(x, \mu, \sigma) \quad F = cdf('norm', x, \mu, \sigma)$
$x = norminv(F, \mu, \sigma) \quad x = icdf('norm', F, \mu, \sigma)$

例9-2 征态分布
征态分布的概率密度函数与分布函数曲线
$(\mu, \sigma^2)参数选择为(-1, 1), (0, 0.1), (0, 1), (0, 10), (1, 1)$

x = [-5: 0.02: 5]';
y1 = [];
y2 = [];
mu1 = [-1, 0, 0, 0, 1];
sig = sqrt([1, 0.1, 1, 10, 1]);
for i = 1: length(mu1)
    y1 = [y1, normpdf(x, mu1(i), sig(i))];
    y2 = [y2, normcdf(x, mu1(i), sig(i))];
end;
plot(x, y1);
figure;
plot(x, y2);

Rayleigh分布
Rayleigh分布的概率密度为:
$p(x) = \left \{ \begin{array}{l} \dfrac{x}{b^2} e^{-\frac{x^2}{2b^2}} \quad x \geq 0 \\ 0 \quad x < 0 \end{array} \right.$
该函数是b的函数
y = raylpdf(x, k), y = pdf(‘rayl’, x, k)
F = raylcdf(x, k), F = cdf(‘rayl’, x, k)
x = raylinv(F, k), x = icdf(‘rayl’, F, k)

例9-7 Rayleigh分布
Rayleigh分布的概率密度函数与分布函数曲线
参数b = 0.5, 1, 3, 5

x = [-eps: -0.02: -0.05, 0: 0.02: 5];
y2 = [];
x = sort(x');
b1 = [.5, 1, 3, 5];
y1 = [];
for i = 1: length(b1)
    y1 = [y1, raylpdf(x, b1(i))];
    y2 = [y2, raylcdf(x, b1(i))];
end;
plot(x, y1);
figure;
plot(x, y2);

随机数与伪随机数生成
为随机数–数学方式生成的，rand(), randn()
生成不同种类分布的随机数的函数调用格式
生成n x m 的 gamma 分布的伪随机数矩阵
A= gamrnd(a, $\lambda$, n, m)
A = random(‘gam’, a, $\lambda$, n, m)
生成$\chi^2$分布的伪随机数矩阵
A = chi2rnd(k, n, m)
A = random(‘chi2’, k, n, m)

09.02 概率计算

离散数据的直方图与饼图表示
一组离散的检测数据 $x_1, x_2, \cdots, x_n$
数据区间(a, b), m等分，$b_1 = a, b_{m+1} = b$
$每个子区间(b_i, b_{i+1}) 落入数据的个数为 f_i, 频度f_i/n$
matlab求解：
k = hist(x, b);
f = k/n;
p = f/($b_2 - b_1$)
直方图: bar(b, f);
饼图:pie(f)

例9-10 Rayleigh分布的直方图
生成30000个伪随机数，满足 Rayleight 分布
参数 b = 1

b = 1;
p = raylrnd(1, 30000, 1);
x = 0: .1: 4;
x1 = x+ 0.05;
yy = hist(p, x1);
yy = yy/(30000*0.1);
bar(x1, yy),
y = raylpdf(x, 1);
line(x, y);
figure;
pie(yy)

例9-11 离散数据的概率表示
随机测得200只荧光灯流明数据，文件c9dlamp.dat,绘图

% A = load('c9dlamp.dat');
A = load('E:\\Work\\201806\\math_4th\\c9dlamp.dat');
bins = [500: 100: 1500] + 50;
f = hist(A, bins)/length(A);
bar(bins, f)
figure;
pie(f)

连续概率问题求解
三个求取概率的公式：
$\xi \leq x$ 的概率：
$P[\xi \leq x] = F(x)$
$x_1 \leq \xi \leq _2$ 的概率：
$P[x_1 \leq \xi \leq x_2] = F(x_2) - F(x_1)$
$\xi \geq x$ 的概率：$P[\xi \geq x] = 1 -F(x)$

例9-12 Rayleigh分布概率计算
已知随机变量x为Rayleigh分布，且b = 1
求出随机变量x值落入区间[0.2, 2]及区间[1, $\infty$)的概率

% 落入区间[0.2, 2] P = F(2) - F(0.2)
b = 1;
p1 = raylcdf(0.2, b);
p2 = raylcdf(2, b);
P1 = p2 - p1

落入区间[1, $\infty$)

p1 = raylcdf(1, b);
P2 = 1 - p1

例9-13 联合概率计算
二维随机变量($\xi, \eta$)的联合概率密度为：
$p(x, y) = \left \{ \begin{array}{l} x^2 + \dfrac{xy}{3}, 0 \leq x \leq 1, 0 \leq y \leq 2 \\ 0, others \end{array} \right.$
求出$P(\xi <1/2, \eta < 1/2) \quad P = \int_{-\infty}^{a} \int_{-\infty}^{b} p(x, y) dy dx$

syms x y;
f = x^2 + x*y/3;
P = int(int(f, x, 0, 1/2), y, 0, 1/2)

例9-14 概率计算
假设某两地A、B间有6个交通岗
在各个交通岗遇到红灯的概率均相同，为p=1/3
中途遇到红灯次数满足二项式分布B(6, p)
试求出从A出发到达B至少遇到一次红灯的概率

x = 0: 6;
y = binopdf(x, 6, 1/3);
P = 1 - y(1) % or P = sum(y(2:end))
p0 = 0.05: 0.05: 0.95;
y = [];
for p = p0
    y = [y 1-binopdf(0, 6, p)];
end;
plot(p0, y, 1/3, P, 'o');

例9-15 利用Monte Carlo 方法计算 $\pi$
试用 Monte Carlo 方法近似求出 $\pi$ 的值 $\pi \approx 4N_1/N$
数学求解公式: $N_1/N \approx \pi / 4$

N = 100000;
x = rand(1, N);
y = rand(1, N);
i = (x.^2 + y.^2) <= 1;
N1 = sum(i);
p = N1/N*4

例9-16 定积分近似
试用 Monte Carlo 法计算积分
$\int_{1}^{3} [1 + e^{-0.2x} \sin (x + 0.5)] dx$
假设 $f(t) \geq 0$
$\dfrac{N_1}{N} \approx \dfrac{1}{M(b-a)} \int_a^b f(x) dx$
$\int_a^b f(x) dx \approx \dfrac{M(b-a) N_1}{N}$
计算公式：$\int_a^b f(x) dx \approx \dfrac{M(b-a)N_1}{N}$

f = @(x)1 + exp(-0.2*x).*sin(x+0.5);
a = 1;
b = 3;
M = 2;
N = 1000000;
x = a+(b-a)*rand(N, 1);
y = M*rand(N, 1);
i = y <= f(x);
N1 = sum(i);
p = M*N1*(b-a)/N,
% 解析解
syms x;
I = vpa(int(1 + exp(-0.2*x)*sin(x+0.5), x, a, b))

例9-17 Brown运动仿真
单个粒子的Brown运动仿真
粒子位置的递推公式
$x_{i+1} = x_i + \sigma \Delta x_i, y_{i+1} = y_i + \sigma \Delta y_i$
运动步距满足标准征态分布
选择比例因子 \sigma = 1

n = 1000;
x = zeros(2, n);
y = zeros(2, n);
s = 1;
r1 = randn(2, n);
for i = 2:n
    x(1, i) = x(1, i-1) + s*r1(1, i);
    y(1, i) = y(1, i-1) +s*r1(2, i);
end;
plot(x(1, :), y(1, :), '-o')

09.03 统计量计算与分析

随机变量的均值与方差
连续随机变量x的概率密度为p(x)
数学期望E[x]:
$E[x] = \int_{-\infty}^{\infty} x p(x) dx$
数学方差D[x] : $D[x] = \int_{-\infty}^{\infty} (x - E[x])^2 p(x) dx$

例9-18 统计量计算
用积分方法求Gamma分布(a > 0, $\lambda$ > 0)的均值与方差
Gamma分布的概率密度 $p(x) = \left \{ \begin{array}{l} \dfrac{\lambda^a x^{a-1}}{\Gamma(a)}e^{-\lambda x}, x \geq 0 \\ 0, \quad x < 0 \end{array} \right.$

syms x;
syms a lam positive;
p = lam^a*x^(a-1)/gamma(a)*exp(-lam*x);
m = int(x*p, x, 0, inf);
s = simplify(int((x-1/lam*a)^2*p, x, 0, inf));
p, m, s
% s需要在早期版本计算，matlab2008a

结果：$m = \dfrac{a}{\lambda}, s = \dfrac{a}{\lambda^2}$

统计量数值计算
在实际中测出一组样本数据 $x = [x_1, x_2, x_3, \cdots, x_n]^T$
它们的均值和方差分别为 – 样本均值方差
$\bar x = \dfrac{1}{n} \sum \limits_{i=1}^{n} x_i, \quad \hat s_x^2 = \dfrac{1}{n} \sum \limits_{i=1}^{n}(x_i-\bar x)^2$
m = mean(x)
s2 = var(x)
无偏方差: $s_x^2 = \dfrac{1}{n-1} \sum \limits_{i=1}^{n}(x_i - \bar x)^2$
s = std(x)
称$s_x \geq 0$ 为 “标准差”

例9-22 伪随机数的统计量
生成一组30000个征态分布随机数
均值为 0.5, 标准差为1.5
分析数据实际的均值，方差和标准差
减小随机变量个数有何影响

p = normrnd(0.5, 1.5, 30000, 1);
mean(p), var(p), std(p)

使用300个随机数：

p = normrnd(0.5, 1.5, 300, 1);
mean(p), var(p), std(p)

随机变量的矩
假设x为连续随机变量, 且p(x)为其概率密度函数
该变量的r阶原点矩
$v_r = \int_{-\infty}^{\infty} x^r p(x) dx$
r阶中心矩
$\mu_r = \int_{-\infty}^{\infty} (x - \mu)^r p(x) dx$
可见, $v_1 = E[x], \mu_2 = D[x]$

例9-21 Gamma分布的矩量计算
考虑Gamma分布(a > 0, \lambda > 0)的原点矩和中心矩

syms x;
syms a lam positive;
p = lam^a*x^(a-1)/gamma(a)*exp(-lam*x);
for n=1:5
    m = int(x^n*p, x, 0, inf)
end;

原点矩通项表达式: $v_k = \dfrac{1}{\lambda^k} \prod \limits_{i=0}^{k-1} (a+i)$

% 直接求出任意阶次的原点矩
syms k;
m = simplify(int((x)^k*p, x, 0, inf))

% 计算中心矩，要用早期版本
for n = 1:7
    s = simplify(int((x-1/lam*a)^n*p, x, 0, inf)),
end;

中心矩没有显而易见的通项公式

矩的数值计算
给定的随机数为一些样本点 $x_1, x_2, \cdots, x_n$
该随机变量的r阶原点矩
$A_r = \dfrac{1}{n} \sum \limits_{i=1}^{n} x_i^r$
matlab: $A_r$ = sum(x.^r)/length(x)

该随机变量的r阶中心矩
$B_r = \dfrac{1}{n} \sum \limits_{i=1}^{n} (x_i - \bar x)^r$
matlab: $B_r$ = moment(x, r)

例9-22 矩量的数值计算
给定一组30000个征态分布随机数
均值为 0.5, 标准差为 1.5, 试求出随机数的各阶矩

A = [];
B = [];
n = 1: 5;
p = normrnd(0.5, 1.5, 30000, 1);
for r = n
    A = [A, sum(p.^r)/length(p)];
    B = [B, moment(p, r)];
end;
A, B

syms x;
A1 = [];
B1 = [];
a = -inf;
b = inf;
p = 1/(sqrt(2*sym(pi))*3/2)*exp(-(x-1/2)^2/(2*(3/2)^2));
for i = 1:5
    A1 = [A1, int(x^i*p, x, a, b)];
    B1 = [B1, int((x-1/2)^i*p, x, a, b)];
end;
A1, B1

09.04 协方差计算

多变量随机数的协方差分析
随机数$(x_1, y_1), (x_2, y_2), (x_3, y_3), \cdots, (x_n, y_n)$
为二维随机变量对(x,y)的样本
二维样本的协方差
$s_{xy} = \dfrac{1}{n-1} \sum \limits_{i=1}^{n} (x_i - \bar x)(y_i - \bar y)$
二维样本的相关系数
$\eta = \dfrac{s_{xy}}{s_x s_y}$

协方差矩阵计算
协方差矩阵
$C = \begin{bmatrix} c_{xx} & c_{xy} \\ c_{yx} & c_{yy} \end{bmatrix}$
其中，$c_{xx} = \sigma_x^2, c_{xy} = c_{yx} = s_{xy}, c_{yy} = \sigma_y^2$
计算协方差矩阵的函数调用格式
C = cov(X)
其中，X的各列均表示不同的随机变量的样本值

例9-23 征态分布的协方差矩阵
使用matlab产生4各满足标准征态分布的随机变量
并求出其协方差矩阵
四列信号相互独立

p = randn(30000, 4);
R = cov(p)

多变量正态分布的联合概率密度及分布函数
给定n组正态分布随机变量 $\xi_1, \xi_2, \cdots, \xi_n$, 它们的均值分别为 $\mu_1, \mu_2, \mu_3, \cdots, \mu_n$, 可以构成一个均值向量 $\mu$, 这些变量的协方差矩阵为 $\Sigma^2$, 这些随机变量的联合概率密度为：
$p(x_1, x_2, \cdots, x_n) = \dfrac{1}{\sqrt{2 \pi}} \Sigma^{-1} e^{-\frac{1}{2}x^T\Sigma^{-2} x}$
其中，$x = [x_1, x_2, \cdots, x_n]^T, \mu = [\mu_1, \mu_2, \cdots, \mu_n]^T$

求随机变量的联合概率密度的函数调用格式
p = mvnpdf(X, $\mu, \Sigma^2$)
p = pdf(‘mvn’, X, $\mu, \Sigma^2$)
其中, X 为n列矩阵
每一列表示一个随机变量

例9-24 协方差矩阵与随机数分布
给定$\mu = [-1, 2]^T, \Sigma^2 = [1, 1; 1, 3]$, 联合概率密度函数

mu1 = [-1, 2];
Sigma2 = [1 1; 1 3];
[X, Y] = meshgrid(-3:0.1:1, -2:0.1:4);
xy = [X(:), Y(:)];
p = mvnpdf(xy, mu1, Sigma2);
P = reshape(p, size(X));
surf(X, Y, P)

若协方差矩阵为对角线矩阵，新概率密度函数

Sigma2 = diag(diag(Sigma2));
p = mvnpdf(xy, mu1, Sigma2);
P = reshape(p, size(X));
surf(X, Y, P)

多变量正态分布随机数生成
产生多变量正态分布随机数的函数调用格式
R = mvnrnd($\mu, \Sigma^2$, m)
R = random(‘mvn’, $\mu, \Sigma^2$, m)
该函数可以生成m组满足多变量正态分布的随机变量，返回的R为m x n 矩阵，每一列表示一个随机变量

例9-24 协方差矩阵与随机数分布
给定$\mu = [-1, 2]^T, \Sigma^2 = [1, 1; 1, 3]$,
二维正态分布的伪随机数的分布情况

mu1 = [-1, 2];
Sigma2 = [1 1; 1 3];
R1 = mvnrnd(mu1, Sigma2, 2000);
plot(R1(:, 1), R1(:, 2), 'o');
figure;
Sigma2 = diag(diag(Sigma2));
R2 = mvnrnd(mu1, Sigma2, 2000);
plot(R2(:, 1), R2(:, 2), 'o');

09.05 离群值检测

离群值、四分位数与盒子图
均值的意义与局限性
例：姚明给一群小朋友讲NBA经历
均值(平均身高)没有意义
中位数(median, 仲数)的定义 $x_1 \leq x_2 \leq \cdots \leq x_n$
n为奇数: $x_{(n+1)/2}$
n为偶数: $(x_{n/2-1} + x_{n/2+1})/2$
matlab求解: median(x)
离群值(outliers):
由中位数将向量x分成两部分
这两部分的中位数称为四分位数(quantile)
第一分位数: $q_1$
第二分位数: $q_2$, 中位数(仲数)
第三分位数: $q_3$
四分位数向量: q = [$q_1, q_2, q_3$]
matlab求解: q = quantile(x, 3)

四分位距(interquartile range, IQR)的定义
第3四分位数与第1四分位数的差: IQR = $q_3 - q_1$
离群值:
正常值的区间: $(q_1 - 1.5 IQR, q_3 + 1.5 IQR)$
正常值区间之外的都是离群值

盒子图绘制：boxplot(x)
上边界与下边界: $(q_1 - 1.5 IQR, q_3 + 1.5 IQR)$

例9-26 盒子图绘制
一组数据，存于c9dlamp.dat文件
四分位数求值及盒子图绘制

% A = load('c9dlamp.dat');
A = load('E:\\Work\\201806\\math_4th\\c9dlamp.dat');
q = quantile(A, 3)
boxplot(A)
sort(A')

离群值检测
单变量数据
Rutgers大学 Niccolo Battistini编写的 outliers()函数
[$v_1, v_2$] = outliers(v, opts, $\alpha$)
选项’grubbs’, quartile’
显著性水平: $\alpha$

function [Y out]=outliers(X,method,alpha)

% Reports vector of observations Y which are not outliers and vector of
% observations out which are outliers according to either Grubbs' test or
% the Quartile method applied to initial vector of observations X
% 
% USAGE:    [Y out]=outliers(X,method,alpha)
% 
% INPUT:
% X         Nx1 vector of observations
% method    string (either 'grubbs' or 'quartile')
% alpha     significance level of Grubbs' test (alpha in [0,1]) or
%           range looseness parameter of Quartile method
% 
% OUTPUT:
% Y         (N-k)x1 vector of observations in X which are not outliers
% out       kx1 vector of observations in X which are outliers

% Author: Niccolo Battistini (Rutgers University)

if nargin~=3;
    error('Wrong # of arguments')
end

switch method

    % Grubbs' test

    case 'grubbs'
        z=1;
        test=false;
        while test==false;
            N=length(X);
            Xmean=repmat(mean(X),N,1); % N by 1 vector of mean of X
            Xstd=std(X); % standard deviation of X
            tcv=tinv((1-alpha)/(2*N),N-2); % t-stat critical value
            Gcv=((N-1)/sqrt(N))*sqrt(tcv^2/(N-2+tcv^2)); % G-stat CV
            [maxdev ind]=max(abs(X-Xmean));
            G=maxdev/Xstd; % G-stat
            if G>Gcv;
                out(z,1)=X(ind);
                X=X([1:N]~=ind);
                z=z+1;
                test=false;
            else
                test=true;
            end
        end
        if z==1;
            out=[];
        end
        Y=X;

    % Quartile method

    case 'quartile'
        Q1=prctile(X,25);
        Q3=prctile(X,76);
        rangedown=Q1-alpha*(Q3-Q1);
        rangeup=Q3+alpha*(Q3-Q1);
        Y=X(X>rangedown & X<rangeup);
        out=X(X<=rangedown | X>=rangeup);
end
end

多变量数据
Antonio Trujillo-Ortiz编写的函数 moutlier1()
moutlier1(X, $\alpha$)

function x = ACR(p,n,alpha);
if nargin < 3,
   alpha = 0.05; %(default)
end; 
if (alpha <= 0 | alpha >= 1)
   fprintf('Warning: significance level must be between 0 and 1\n');
   return;
end;
if nargin < 2, 
   error('Requires at least two input arguments.');
   return;
end;
a = alpha;
Fc = finv(1-a/n,p,n-p-1); %F distribution critical value with p and n-p-1 degrees of freedom using the Bonferroni correction. 
ACR = (p*(n-1)^2*Fc)/(n*(n-p-1)+(n*p*Fc)); % = ((-1*((1/(1+(Fc*p/(n-p-1))))-1))*((n-1)^2))/n;
x = ACR;
return,

function moutlier1(X,alpha)
if nargin < 2, 
   alpha = 0.05;  %(default)
end 

if nargin < 1, 
   error('Requires at least one input arguments.');
end

mX = mean(X); %Means vector from data matrix X.
[n,p] = size(X);
difT = [];

for j = 1:p;
   eval(['difT=[difT,(X(:,j)-mean(X(:,j)))];']); %squared Mahalanobis distances.
end

S = cov(X);
D2T = difT*inv(S)*difT'; 
[D2,cc] = sort(diag(D2T));  %Ascending squared Mahalanobis distances.

D2C = ACR(p,n,alpha);

idx = find(D2 >= D2C);
o = cc(idx);
io = D2(idx);

if isempty(o);
    disp(' ')
    fprintf('With a given significance level of: %.2f\n', alpha);
    disp('Non observation(s) resulting as multivariate outlier(s).');
else
    disp(' ')
    disp('Table of observation(s) resulting as multivariate outlier(s).')
    fprintf('----------------------------------------------\n');  
    disp('                            D2');
    disp('Observation              observed');
    fprintf('----------------------------------------------\n');  
    fprintf(' %6.0f               %10.4f\n',[o,io].');
    fprintf('----------------------------------------------\n');  
    fprintf('With a given significance level of: %.2f\n', alpha);
    fprintf('Critical value for the maximum squared Mahalanobis distance: %.4f\n', D2C);
    disp('D2 = squared Mahalanobis distance.');
end

return,

例9-27 单变量离群值
数据表格，存于c9dlamp.dat文件
离群值检测

% A = load('c9dlamp.dat');
A = load('E:\\Work\\201806\\math_4th\\c9dlamp.dat');
[v1 v2] = outliers(A, 'grubbs', 0.05) % 95%的显著性水平

例9-28 多变量离群值检测
NBA球队数据

% 数据读入
X = [447, 149, 22.8; 401, 160, 13.5; 356, 119, 49; 338, 117, -17.7; 328, 109, 2;
    290, 97, 25.6; 284, 102, 23.5; 283, 105, 18.5; 282, 109, 21.5; 280, 94, 10.1;
    278, 82, 15.2; 275, 102, -16.8; 274, 98, 28.5; 272, 97, -85.1; 258, 72, 3.8; 
    249, 96, 10.6; 244, 94, -1.6; 239, 85, 13.8; 236, 91, 7.9; 230, 85, 6.9;
    227, 63, -19.7; 218, 75, 7.9; 216, 80, 21.9; 208, 72, 15.9; 202, 78, -8.4;
    199, 80, 13.1; 196, 70, 2.4; 188, 70, 7.8; 174, 70, -15.1];
% 离群值检测
moutlier1(X, 0.05)

检测结果图形表示
由检测结果可知，第14队数据是离群值

% 三维散点图显示
plot3(X(:, 1), X(:, 2), X(:, 3), 'o');
% x-z平面图
figure;
plot(X(:, 1), X(:, 3), 'o');
% y-z 平面图
figure;
plot(X(:, 2), X(:, 3), 'o');

09.06 参数估计与区间估计

参数估计与区间估计
求取参数与区间估计的函数调用格式
[$\mu, \sigma^2, \Delta \mu, \Delta \sigma^2$] = normfit(x, $P_{ci}$)
[$\mu, \sigma^2, \Delta \mu, \Delta sigma^2$] = fittest(‘norm’, x, $P_{ci}$)
其中，$x = [x_1, x_2, \cdots, x_n]^T$ 是实测一组数据
$\mu 是该分布的均值, \sigma^2 是该分布的方差$
$\Delta \mu 及 \Delta \sigma^2 是置信区间$
$P_{ci} 为用户指定的置信度$

给定分布的均值与方差计算
函数norminv()可用于求出相关值，这样就可以得出所需的参数
Gamma分布的参数($a, \lambda$): gamfit()
Rayleigh分布的参数估计函数为: raylfit()
均与分布的参数估计函数为: unifit()
Poisson分布的参数估计函数为: poissfit()
可以调用 fittest 函数

例9-29 Gamma分布的参数估计
试用 gamrnd() 函数生成一组 a = 1.5, $\lambda$ = 3 的伪随机数
用参数估计的方法以不同的置信度进行估计
选择置信度为 90%, 92%, 95%, 98%
置信度，估计值，估计区间

p = gamrnd(1.5, 3, 30000, 1);
Pv = [0.9, 0.92, 0.95, 0.98];
A = [];
for i = 1: length(Pv)
    [a, b] = gamfit(p, Pv(i));
    A = [A; Pv(i), a(1), b(:,1)', a(2), b(:, 2)'];
end;
A

多元线性回归与区间估计
输出信号y
n路输入信号: $x_1, x_2, \cdots, x_n$
线性组合
$y = a_1x_2 + a_2x_2 + a_3x_3 + \cdots + a_nx_n$
其中， $a_1, a_2, \cdots, a_n$ 为待定系数
线性回归
求解线性代数方程

得到的实测数据满足的关系式
实测数据
$y_1 = x_{11}a_1 +x_{12}a_2 + \cdots + x_{1n}a_n + \varepsilon_1$
$y_2 = x_{21}a_2 +x_{22}a_2 + \cdots + x_{2n}a_n + \varepsilon_2$
$\vdots$
$y_m = x_{m1}a_1 +x_{m2}a_2 + \cdots + x_{mn}a_n + \varepsilon_m$
观测数据组与未知待定参数个数不同
尽量使得总体误差最小

参数估计的最小二乘求解
目标函数选择为使得残差的平方和最小:
$J = \min \limits_{a} \varepsilon^T \varepsilon$
系数向量a为 $\hat a = (X^T X)^{-1} X^T y$
求最小二乘解的函数调用格式
$a = X \ y \quad a = inv(X'*X)*X'*y$
求解函数, 1-a为用户指定的置信度
[$\hat a, a_{ci}$] = regress(y, X, a)

例9-30 参数估计与区间估计
给定线性回归方程如下，生成120组随机输入值$x_i$
计算输出向量y，估计出系数$a_i$
$y = x_1 - 1.232 x_2 + 2.23 x_3 + 2x_4 + 4 x_5 + 3.792x_6$
用最小二乘计算公式

a = [1, -1.232, 2.23, 2, 4, 3.792]';
X = randn(120, 6);
y = X*a;
a1 = inv(X'*X)*X'*y
% 计算出98%的置信度区间
[a, aint] = regress(y, X, 0.02)

混叠噪声的信号参数估计
给输出样本叠加N(0, 0.5)区间的正态分布噪声
再绘制参数估计的置信区间

yhat = y + sqrt(0.5)*randn(120, 1);
[a, aint] = regress(yhat, X, 0.02)
errorbar(1:6, a, aint(:,1)-a, aint(:,2)-a)

将噪声方差设为0.1

yhat = y + sqrt(0.1)*randn(120, 1);
[a, aint] = regress(yhat, X, 0.02)
errorbar(1:6, a, aint(:, 1)-a, aint(:, 2)-a)

非线性函数的最小二乘参数估计与区间估计
假设数据 $x_i, y_i, i = 1, 2, \cdots, N$满足原型函数
$\hat y(x) = f(a, x)$
原函数严格写成
$\hat y(x) = f(a, x) + \varepsilon$
引入目标函数:
$I = \min \limits_{a} \sum \limits_{i=1}^{N} [y_i - \hat y(x_i)]^2 = \min \limits_{a} \sum \limits_{i=1}^{N}[y_i - f(a, x_i)]^2$

参数估计的函数调用格式
最小二乘拟合： [a, r, J] = nlinfit(x, y, fun, a0)
由置信度为95%的置信区间 ： c = nlparci(a, r, J)
与函数lsqcurvefit()的功能相似

例9-31 参数与区间估计
原型函数 $y(x) = a_1 e^{-a_2x} + a_3e^{-a_4x} \sin(a_5x)$
95%置信度的置信区间
叠加均匀分布的噪声信号再进行参数与区间估计

f = @(a, x)a(1)*exp(-a(2)*x) + ...
    a(3)*exp(-a(4)*x).*sin(a(5)*x);
x = 0: 0.1: 10;
y = f([0.12, 0.213, 0.54, 0.17, 1.23], x);
[a, r, j] = nlinfit(x, y, f, [1;1;1;1;1]);
a
ci = nlparci(a, r, j)

叠加噪声后估计
样本点数据 $y_i$
叠加上[0, 0.02] 区间均匀分布的噪声信号

y = f([0.12, 0.213, 0.54, 0.17, 1.23], x) ...
    + 0.02*rand(size(x));
[a, r, j] = nlinfit(x, y, f, [1;1;1;1;1]),
ci = nlparci(a, r, j)
errorbar(1:5, a, ci(:, 1)-a, ci(:, 2)-a)

09.07 统计假设检验

统计假设检验的概念及步骤
为什么要假设检验？
生产一批灯泡，标明寿命，怎么标？
产品设计指标(预期寿命)
随机选n个产品，测出相关的参数
检验是否合格
假设检验的一般步骤
先假设总体具有某种统计特征
然后再检验这个假设是否可信
这种方法称为“统计假设检验方法”
统计假设检验在统计学中是有重要地位的

假设检验的步骤
做出假设–产品合格: $\mathscr{H_0}: \mu = \mu_0$
随机选择n个产品，测出均值$\bar x$ 与标准差s
构造统计量 $\mu = \sqrt{n}(\bar x - \mu_0)/s$, 满足N(0, 1)分布
求出逆概率分布$\int_{K_{\alpha/2}}^{K_{\alpha/2}} \dfrac{1}{\sqrt{2 \pi}} e^{x^2/2} dx < 1 - \alpha$
$K_{\alpha/2}$ = norminv(1 - $\alpha$/2, 0, 1)
$K_{\alpha/2}$ = icdf(‘norm’, 1 - $\alpha/2$, 0, 1)
做出结论：若|$\mu$| < $K_{\alpha/2}$, 不能拒绝假设

例9-34 工艺变化的强度检验
已知某产品的平均强度 $\mu_0 = 9.94kg$
现改变工艺，在新产品中随机抽取200件
平均强度为: $\bar x = 9.73kg$
标准差为: $s = 1.62kg$
问改变工艺对产品强度有无显著影响
引入两个命题:
$\left \{ \begin{array}{l} \mathscr{H}_0 : \mu = \mu_0 \; no significant change \\ \mathscr{H}_1: reject \; \mathscr{H}_0 \end{array} \right.$

选取统计量 $\mu = \dfrac{\sqrt{n}(\bar x - \mu_0)}{s}$
$\int_{K_{\alpha/2}}^{K_{\alpha/2}} \dfrac{1}{\sqrt{2 \pi}} e^{x^2/2} dx < 1 - \alpha \quad |\mu| < K_{\alpha/2}$
该统计量满足标准正态分布N(0, 1)

n = 200;
mu0 = 9.94;
xbar = 9.73;
s = 1.62;
u = sqrt(n)*(mu0 - xbar)/s
alpha = 0.02;
K = norminv(1 - alpha/2, 0, 1),
H = abs(u) < K

两组数据是否有明显差异
有两组数据，第一组随机选$n_1$个样本，第二组$n_2$个样本
假设检验的步骤
做出假设(没有差异) – $\mathscr{H_0}: \mu_1 = \mu_2$
构造T分布统计量 $t = \dfrac{\bar x_1 - \bar x_2}{\sqrt{s_1^2/n_1 + s_2^2/n_2}}$
k = min($n_1 - 1, n_2 - 1$)
$T_0 = \mathbb{tinv}(\alpha/2, k), or \; T_0 = \mathbb{icdf}('t', \alpha/2, k)$
若|t| < |$T_0$|, 不能拒绝

例9-35 失眠药物的药效评价
失眠病患者，随机分成两组， 各10人
不同组使用不同药物，测延长睡眠小时数
$\begin{array}{l| l l }\hline A & 1.9 & 0.8 & 1.1 & 0.1 & -0.1 & 4.4 & 5.5 & 1.6 & 4.6 & 3.4 \\ B &0.7 & -1.6 & -0.2 & -1.2 & -0.1 & 3.4 & 3.7 & 0.8 & 0 & 2 \\ \hline \end{array}$
两种药物是否有显著差异？

x = [1.9, 0.8, 1.1, 0.1, -0.1, 4.4, 5.5, 1.6, 4.6, 3.4];
y = [0.7, -1.6, -0.2, -1.2, -0.1, 3.4, 3.7, 0.8, 0, 2];
n1 = length(x);
n2 = length(y);
k = min(n1-1, n2-1);
t = (mean(x) - mean(y)) ...
    /sqrt(std(x)^2/n1 + std(y)^2/n2)
a = 0.05;
T0 = tinv(a/2, k),
H = abs(t) < abs(T0)
% 图解分析
boxplot([x.', y.'])

正态分布的均值假设检验
已知一组数据
该数据复合正态分布，且已知其标准差为$\sigma$
假设其均值为$\mu$
如何检验？$\mathscr{H_0}: \mu = \mu_0$
[H, s, $\mu_{ci}$] = ztest(X, $\mu_0, \sigma, \alpha$)
如果未知数据的标准差，如何检验？
[H, s, $\mu_{ci}$] = ttest(X, $\mu_0, \alpha$)

例9-36 正态分布的均值假设检验
生成一组400个正态分布随机数
均值为1，标准差为2
由数据检验以下若标准差$\sigma$-2,数据的均值是否为1

r = normrnd(1, 2, 400, 1);
[H, p, ci] = ztest(r, 1, 2, 0.02)
% 如果未知标准差，检验均值是否为1
[H, p, ci] = ttest(r, 1, 0.02)

正态性假设检验
测得一组数据，检验其是否满足正态分布
Jarque-Bera假设检验
[H, s] = jbtest(X, $\alpha$)
Lilliefors假设检验
[H, s] = lillietest(X, $\alpha$)
如果满足正态分布，则找出均值与方差
[$\mu, \sigma, \mu_{ci}, \sigma_{ci}$] = normfit(X, $\alpha$)
图形验证： normplot(X)

例9-37 判定分布是否为正态分布
前面给出的表格，文件 c9dlumen.dat
正态性判定，并得出正态分布参数

% X = load('c9dlumen.dat');
X = load('E:\\Work\\201806\\math_4th\\c9dlumen.dat');
[H, p] = jbtest(X, 0.05)
[mu1, sig1, mu_ci, sig_ci] = normfit(X, 0.05)
normplot(X)

例9-38 非正态分布的检验结果
生成一组400个Rayleigh分布数据
正态性检验结果

X = raylrnd(1.5, 400, 1);
[H, p, c, d] = jbtest(X, 0.05)
% 图形验证
normplot(X)

任意分布的Kolmogorov-Smirnov检验
Kolmogorov-Smirnov检验是检验任意已知分布函数的一种有效的假设检验算法
函数调用格式: [H, s] = kstest(X, cdffun, $\alpha$)
其中，cdffun为两列均值, 第一列为自变量，第二列应该为要检验的分布函数在自变量处的值
某些分布的图形验证: probplot()

例9-40 Rayleigh分布的检验
生成一组Rayleigh分布数据
对生成的随机数进行假设检验，是否满足分布
生成Rayleigh分布的数据

r = raylrnd(1.5, 400, 1);
b = raylfit(r)
% 假设检验
s = sort(r);
[H, p] = kstest(s, [s raylcdf(s, b)], 0.05)
% 图形验证
probplot('rayleigh', r)

09.08 方差分析

方差分析
方差分析(analyysis of variance, ANOVA)是应该遗传学家，统计学家 Ronald Fischer 提出的一种分析方法
方差分析技术是假设检验的拓展
$\mathscr{H_0} : \mu_1 = \mu_2 = \cdots = \mu_N$
单因子方差分析
双因子方差分析
多因子方差分析

单因子方差分析
单因子方差分析就是指对一些观察来说，只有一个外界因素可能对观测的现象产生影响
求解单因子方差分析的函数调用格式:
[p, tab, stats] = anova1(X)
其中，X为需要分析的数据
若 p < $\alpha$, 则拒绝假设
$\mathscr{H_0}: \mu_1 = \mu_2 = \cdots = \mu_N$

单因子方差分析表
ANOVA表
$\begin{array}{|l| l } \hline \mathbb{source} & \mathbb{sum \; of \; squares} & DOF & \mathbb{mean \; squares} & F & \mathbb{probability \; p} \\ \hline groups & SSA = \sum \limits_{i} n_i\bar y_{i,:}^2 - N \bar y_{:,:}^2 & I-1 & MSSA = \dfrac{SSA}{I-1} & \dfrac{MSSA}{MSSE} & p = P(F_{I-1, N-I} > c) \\ \hline error & SSE = \sum \limits_{i} \sum \limits_{k} y_{i,k}^2 - \sum \limits_{i} n_i \bar y_{i,:}^2 & N-I & MSSE = \dfrac{SSE}{N-I} & & \\ \hline total & SST = \sum \limits_{i} \sum \limits_{k} y_{i,k}^2 - N \bar y_{:,:}^2 & N-1 & & & \\ \hline \end{array}$
p的值是不是很小
盒子图的观察

例9-42 单元素方差分析
有5种药物比较疗效
将30个病人随机地分成5组
每组使用同一种药物，并记录病人治疗时间
评价疗效–5种药物疗效是否有显著不同?

A = [5, 4, 6,7,9; 8, 6, 4, 4, 3; ...
    7,6,4, 6, 5; 7, 3, 5, 6, 7; ...
    10, 5, 4, 3, 7; 8, 6, 3, 5, 6];
mean(A),
[p, tbl, stats] = anova1(A)

盒子图的直接观察

双因子方差分析
如果有两种因子可能影响到某现象的统计规律，则应该引入双因子方差分析的概念
观测量y可以表示为一个三维数组 $y_{i,j,k}$, 表示第1个因子取第i个水平，第二个因子取第j个水平时，组内第k个对象的观测指标
三个假设：
$\alpha_i为第一因子单独作用对现象没有影响$
$\mathscr{H_1}: \alpha_1 = \alpha_2 = \cdots = \alpha_I$
$\beta_j 为第二因子单独作用对现象没有影响$
$\mathscr{H_2}: \beta_1 = \beta_2 = \cdots = \beta_J$
$\gamma_k 为两个因子同时作用效应$
$\mathscr{H_3}: \gamma_1 = \gamma_2 = \cdots = \gamma_{IJ}$
三个概率的定义及意义
$若 p_A < c_1 则拒绝假设 \mathscr{H_1}$
$p_A = P(F_{[I-1, IJ(K-1)]} > c_1)$
$若p_B < c_2, 则拒绝假设 \mathscr{H_2}$
$p_B = P(F_{[J-1, IJ(K-1)]} > c_2)$
$若p_{AB} < c_3, 则拒绝假设 \mathscr{H_3}$
$p_{AB} = P(F_{[(I-1)(J-1), IJ(K-1)]} > c_3)$
双因子方差分析：[p, tab, stats] = anova2(X)

双因子方差表
$\begin{array}{l|l}\hline source & square \; of \; sums & DOF & mean \; squared \; error & F & p \\ \hline factor \; A & SSA = JK\sum \limits_{i} \bar y_{i,:,:}^2 - IJK \bar y_{:,:,:}^2 & I-1 & MSSA = \dfrac{SSA}{I-1} & \dfrac{MSSA}{MSSE} & p_A \\ factor\;B & SSB = IK\sum \limits_{i} \bar y_{:,j,:}^2 - IJK \bar y_{:,:,:}^2 & J-1 & MSSB = \dfrac{SSB}{J-1} & \dfrac{MSSB}{MSSE} & p_B \\ interation & SSAB = K \sum \limits_{ij} \bar y_{i,j,:}^2 - JK \sum \limits_{i} \bar y_{i,:,:}^2 + IJK \bar y_{:,:,:}^2 & (I-1)(J-1) & MSSAB = \dfrac{SSAB}{(I-1)(J-1)} & \dfrac{MSSAB}{MSSE} & p_{AB} \\ errors & SSE = \sum \limits_{ijk} y_{i,j,k}^2 - K\sum \limits_{i} \sum \limits_{j} \bar y_{i, j,:}^2 & IJ(K-1) & MSSE = \dfrac{SSE}{IJ(K-1)} & & \\ totoal & SST = \sum \limits_{ijk} y_{i,j,k}^2 - IJK \bar y_{:,:,:}^2 & IJK-1 & & & \\ \hline \end{array}$

例9-43 树的生长地、树种的影响
比较3种松树在4个不同地区的生长情况有无差别
在每个地区对每种松树随机地选择5株
测量它们的胸径

B = [23, 15, 26, 13, 21, 25, 20, 21, 16, 18, 21, 17, 16, 24, 27, 14, 17, 19, 20, 24;
    28, 22, 25, 19, 26, 30, 26, 26, 20, 28, 19, 24, 19, 25, 29, 17, 21, 18, 26, 23;
    18, 10, 12, 22, 13, 15, 21, 22, 14, 12, 23, 25, 19, 13, 22, 16, 12, 23, 22, 19];
anova2(B', 5);
[p a b] = anova2(B', 5);
% 均值计算
C = [];
for i = 1: 3
    for j = 1: 4
        C(i, j) = mean(B(i, [1:5] + (j-1)*5));
    end;
end;
C = [C; mean(C)];
C = [C mean(C')']

多因子方差分析
manoval()函数进行多因子方差分析

09.09 主成分分析

主成分分析
假设某一事件发生可能受$x_1, x_2, \cdots, x_N$因素影响，识别出到底哪些因素起主要作用，哪些可以忽略
而实测数据共有M组
这样可以假设这些数据由一个 $N \times M$ 矩阵X表示
记该矩阵的每一列的均值为 $\bar x_i, i = 1, 2, \cdots, N$
mean(X)
由X建立起协方差矩阵R
$r_{ij} = \dfrac{\sqrt{\sum \limits_{k=1}^{M}(x_{ki} - \bar x_i)(k_{kj} - \bar x_j)}}{\sqrt{\sum \limits_{k=1}^{M}(x_{ki} - \bar x_{i})^2 \sum \limits_{k=1}^{M}(x_{kj} - \bar x_j)^2}}$
R = corr(X -mean(X))
由R矩阵可以分别得出
特征向量$e_i$
对应的排序特征值 $\lambda_i, \lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_N \geq 0$
特征向量矩阵的每一列进行相应的归一化
$\Vert e_i \Vert = 1 或 \sum \limits_{j=1}^{N} e_{ij}^2 = 1$
[e, d] = eig(R); d = diag(d);
d = d(end: -1: 1); e = fliplr(e);

找出哪些因素贡献大
计算主成分贡献率和累计贡献率
主成分贡献率: $\gamma_i = \dfrac{\gamma_i}{\sum \limits_{k=1}^N \lambda_k}$
累计贡献率: $\delta_i = \dfrac{\sum \limits_{k=1}^{i} \lambda_k}{\sum \limits_{k=1}^N \lambda_k}$
若前n个特征值的累计贡献率大于 85% ~ 95%
可以认为这n个因素是原问题的主成分
d/sum(d)

构造变换矩阵
建立新变量坐标 Z = XL,即
$\left \{ \begin{array}{l} z_1 = l_{11}x_1 + l_{21}x_2 + \cdots + l_{N1}x_N \\ z_2 = l_{12}x_1 + l_{22}x_2 + \cdots + l_{N2}x_N \\ \vdots \\ z_N = l_{1N}x_1 + l_{2N}x_2 + \cdots + l_{NN}x_N \end{array} \right.$
其中变换矩阵第i列的系数 $l_{ji}$ 可以计算 $l_{ji} = \sqrt{\lambda_i} e_{ji}$
D = [d’; d’; \cdots; d’];
L = real(sqrt(D)).*e, Z = X*L

降维矩阵的构造
若前n个成分做主成分，则矩阵L的m列以后各值应该趋于0，上式化为
$\left \{ \begin{array}{l} z_1 = l_{11}x_1 + l_{21}x_2 + \cdots + l_{N1}x_N \\ \cdots \\ z_n = l_{1n}x_1 + l_{2n}x_2 + \cdots + l_{Nn}x_N \end{array} \right.$
即，在适当的线性变换下，原来的N维问题就可以简化成n维问题

例9-44 主成分分析的降维处理
假设某三维曲线上的样本点由下列函数直接生成
$x = t \cos 2t, y = t \sin 2t, z = 0.2x + 0.6y$
试用主成分分析的方法对其降维处理

t = [0: 0.1: 3*pi]';
x = t.*cos(2*t);
y = t.*sin(2*t);
z = 0.2*x + 0.6*y;
X = [x y z];
R = corr(X);
[e, d] = eig(R),
d = diag(d),
plot3(x, y, z)

% 降维处理
d = d(end:-1:1);
e = fliplr(e);
D = [d'; d'; d'];
L = real(sqrt(D)).*e,
Z = X*L;
% 降维效果
plot(Z(:, 1), Z(:,2))

降维矩阵与坐标变换(Z =XL)
$L=\begin{bmatrix}0.184 & -0.9829 & 0 \\ 0.9653 & 0.261 & 0 \\ 0.9975 & -0.0713 & 0 \end{bmatrix} \Rightarrow \left \{ \begin{array}{l} z_1 = 0.1840 x + 0.9653y + 0.9975z \\ z_2 = -0.9829 x + 0.2610 y - 0.0713 z \end{array} \right.$