matlab/Matlab概率统计与曲线拟合

一、二项分布

  • 二项分布来源于伯努利试验 (事件发生概率 $\boldsymbol{p}$ ) :

​ 含义为独立重复N次试验后, 事件总共发生k次的概率

  • 分布函数 $F(X=k)=P\{X \leq k\},$ 二项分布记为 $B(n, p)$
  • $\mathrm{p_k}=$ binopdf $(\mathrm{k}, \mathrm{N}, \mathrm{p}) \quad$ 获得事件共发生$k$次的概率 $\boldsymbol{P}\{\boldsymbol{x}=\boldsymbol{k}\}$
  • $\mathrm{p_k}=$ binocdf $(\mathrm{k}, \mathrm{N}, \mathrm{p}) \quad$ 为事件最多发生$k$次的概率 $\boldsymbol{P}\{\boldsymbol{x} \leq \boldsymbol{k}\}$
  • $\mathrm{R}=$ binornd $(\mathrm{N}, \mathrm{p}, \mathrm{m}, \mathrm{n}) \quad$ 将生成一个服从二项分布 $B(N, \boldsymbol{p}),$
    规模为 $m \times n$ 的随机矩阵
  • 二项分布的数字特征 $E(k)=N p, D(k)=N p(1-p)$

例:画出$N=100,p=0.5$情况下的二项分布概率特性曲线

1
2
3
4
5
N = 100; p = 0.5;					% 总试验次数和单次试验发生概率
k = 0:N; % 所有可能的事件发生次数
pdf = binopdf(k, N, p); % 绘制概率曲线
cdf = binocdf(k, N, p); % 绘制分布曲线
h = plotyy(k, pdf, k, cdf); % 左右两侧不同的纵轴刻度代表两个函数

绘制结果为:

image-20201124102008010

进阶绘图技巧:set函数的使用

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
set(get(h(1),'Children'),'Color','b','Marker', '.', 'MarkerSize', 13)
% 句中 get(h(1),'children') 表示获取刚才第一条曲线绘制的所有子对象
% 然后将第一条曲线改为蓝色 并且在采样点加注实心点 不是只画点

set(get(h(1),'Ylabel'), 'String','pdf')
% 句柄包含多个绘图时 需要 get 出来再操作 此行改变了左侧 Y 轴的标记名

set(h(2),'Ycolor',[1, 0, 0])
% 第二条曲线纵坐标轴颜色改为纯红色

set(get(h(2),'Children'),'Color','r','Marker','+','MarkerSize',4)
% 第二条曲线改为红色并用来标注采样点

set(get(h(2), 'Ylabel'),'String', 'cdf')
% 右侧Y轴的标记名
xlabel('k')
grid on

绘图结果为:

image-20201124103134184

二、正态分布

  • 正态分布 $N\left(\mu, \sigma^{2}\right)$ 为连续型随机分布,期望 $\mu,$ 标准差 $\sigma$ :

对应分布函数 $\Phi(x)=\int_{-\infty}^{x} f(t) \mathrm{d} t$

  • $\mathrm{px}=$ normpdf $(\mathrm{x}, \mu,$ $\sigma$) $\quad$ 获得服从正态分布 $N(\mu, \sigma^2)$ 的随机变量概率密度函数在$x$的取值。
  • $\mathrm{px}=$ normcdf $(\mathrm{x}, \mu,$ $\sigma$) $\quad$ 获得上述正态分布随机变量不超过$x$的总概率

  • $\mathrm{R}=$ normrnd $(\mu, $ $\sigma$ $, \mathrm{m}, \mathrm{n})$ 将生成一个服从上述正态分 布,规模为 $m \times n$ 的随机样本构成的矩阵

  • $\mathrm{R}=$ randn $(\mathrm{m}, \mathrm{n})$ 将生成一个服从标准正态分布 $N(\mathbf{0}, \mathbf{1}),$ 规模为 $m \times n$ 的随机样本构成的矩阵,事实上,我们可以利用这个矩阵可以构造任何正态分布随机矩阵。

例:正态分布几何表示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
mu = 3; sigma = 0.5;
x = mu + sigma*[-3:-1,1:3]; % 设置六个不同的采样点
yf = normcdf(x, mu, sigma); % 获得六个点的cdf值
P = [yf(4)-yf(3), yf(5)-yf(2), yf(6)-yf(1)]; % 计算cdf的差值(内部区域面积)
xd = 1:0.1:5; yd = normpdf(xd, mu, sigma); clf

for k=1:3
xx = x(4-k):sigma/10:x(3+k);
yy = normpdf(xx, mu, sigma);
% 对于三个不同的面积区间进行不同范围的采样,并获得 pdf 函数的值
subplot(3, 1, k), plot(xd, yd, 'b'); % 绘图位于3行1列第k个位置
hold on, fill([x(4-k), xx, x(3+k)], [0, yy, 0], 'g'); hold off
if k<2
text(3.8, 0.6, '[{\mu}-{\sigma}, {\mu}+{\sigma}]')
else
kk = int2str(k);
text(3.8, 0.6, ['[{\mu}-', kk, '{\sigma}, {\mu}+', kk, '{\sigma}]'])
end
text(2.8, 0.3, num2str(P(k))); shg % 填充区域内显示面积
end
xlabel('x'); shg

绘图结果为:

image-20201124105059445

三、统计分析命令

  • $\min (\mathrm{x}), \max (\mathrm{x})$ 分别计算矩阵各列的最大值或最小值,
  • 若计算整个矩阵最大元素, 可用$\max (\max (\mathrm{x}))$ 或$\max (\mathrm{X}(:))$
  • $\mathrm{mean} (\mathrm{x}),$ $\mathrm{median} (\mathrm{x})$ 分别计算矩阵各列的均值与中位值
  • $\mathrm{std}(\mathrm{x})$,$ \mathrm{var} $$(\mathrm{x})$ 分别计算矩阵各列的样本标准差与样本方差
  • $\mathrm{cov}$ $(\mathrm{x})$ 计算矩阵$x$各列所组成列向量计算出的协方差矩阵, 注意到对应的分丹仍然是$m-1$
  • $\mathrm{corrcoef} (\mathrm{x}) $计算矩阵x各列所组成列向量对应的相关系数

例:产生1000个服从$N(2, 0.5^2)$的随机数

1
2
3
4
5
6
7
8
mu = 2; s = 0.5;
rng(22, 'v5normal')
x = randn(1000, 1);

y = s*x+mu;

subplot(2, 1, 1), histfit(x), axis([-5, 5, 0, 100]), ylabel('x')
subplot(2, 1, 2), histfit(y), axis([-5, 5, 0, 100]), ylabel('y')

4

四、多项式拟合

假设 $y=a_{1} x^{n}+a_{2} x^{n-1}+\cdots+a_{n} x+a_{n+1},$ 我们获取其函数曲线上的一组采样点 $\left\{x_{i}, y_{i} \mid i=1,2, \ldots, m\right\},$ 利用数学方法确定或估计系数 $a_{1}, a_{2}, \ldots, a_{n+1}$ 的问题称之为多项式拟合问题

一般来讲,多项式拟合往往会与逼近或插值这两种知识相结合。在采样点准确,函数光滑的情况下,高阶的拟合(即假设更大的 $n$ ) 往往效果更佳,但如果采样信息有噪声误差, 过大的$n$可能会让结果失去拟合意义(一般设定 $n<m$)

  • p = ployfit(x, y, n)将通过数组$x$和$y$的数据进行拟合,拟合的阶数或次数设定为自然数$n$,返回多项式系数$p$
  • yy=polyval (p, x)可以将多项式系数回代,观察拟合值

利用MATLAB函数计算采样$y$值向量与拟合值向量误差的$2-$范数, $1-$范数与$\infty-$范数,可以分析其平方残差绝对值残差一致逼近残差的情况。

例:多项式拟合实例

1
2
3
4
5
6
7
x0 = 0:0.1:1;
y0 = [-.447, 1.978, 3.11, 5.25, 5.02, 4.66, 4.01, 4.58, 3.45, 5.35, 9.22]; % 构造原始数据
n = 3;P = polyfit(x0, y0, n) % 多项式拟合
xx = 0:0.01:1; yy = polyval(P, xx); % 利用得到的多项式代回得到预测值
plot(xx, yy, '-b', x0, y0, '.r', 'MarkerSize', 20); % 绘图
legend('拟合曲线', '原始数据', 'Location', 'SouthEast')
xlabel('x')

绘制结果为:

5

进阶制表

1
2
y1 = polyval(P, x0);
T = table(x0', y0', y1', y1'-y0', 'VariableNames', {'X', 'Y', 'Fit', 'FitError'})

表打印结果为:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
T =

11×4 table

X Y Fit FitError
___ ______ ________ _________

0 -0.447 -0.90431 -0.45731
0.1 1.978 2.2819 0.3039
0.2 3.11 4.0659 0.95592
0.3 5.25 4.7879 -0.46211
0.4 5.02 4.788 -0.23204
0.5 4.66 4.4063 -0.25372
0.6 4.01 3.983 -0.027002
0.7 4.58 3.8583 -0.72174
0.8 3.45 4.3722 0.92223
0.9 5.35 5.865 0.51503
1 9.22 8.6768 -0.54316

1. 多项式拟合的最小二乘理解

polyfit函数的方法即解最小二乘问题 :

方法是构造 $y=\left[y_{1}, y_{2}, \ldots, y_{m}\right]^{\top}, a=\left[a_{1}, a_{2}, \ldots, a_{n+1}\right]^{\top}$$,X=\left[\begin{array}{cccc}x_{1}^{n} & \cdots & x_{1} & 1 \\ x_{2}^{n} & \cdots & x_{2} & 1 \\ \vdots & \vdots & \vdots & \vdots \\ x_{m}^{n} & \cdots & x_{m} & 1\end{array}\right],$

易得 $y=X a,$ 在 $m>n+1$ 时, 方程超定 (可能无解) , 此时最小二乘解可以通过 $a=x\backslash y$ 获得

例:用最小二乘法获得拟合结果

1
2
3
4
5
6
7
8
9
10
11
12
13
x0 = (0:0.1:1)';
y0 = [-.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22]';

m = length(x0);
n = 3;
X = zeros(m,n+1); %m个采样点,n+1个未知系数

for k=1:n
X(:, n-k+1) = (x0.^k);
end

X(:, n+1) = ones(m, 1);
aT = (X\y0)'

输出结果应与调用多项式拟合函数得到的P相同

1
2
3
aT =

56.6915 -87.1174 40.0070 -0.9043

2. 适用于特定问题的拟合或回归方法

  • 岭回归模型:本质上仍可化归为最小二乘问题

其中, $\boldsymbol{a}=\left[\begin{array}{llll}a_{1} & a_{2} & \cdots & a_{n}\end{array}\right]^{\top}$ 表示拟合系数。

  • Lasso模型:对拟合系数具有稀疏正则的模型:
  • 最小绝对残差 (LAR)模型:对离群值的处理有更好效果 :

实际问题中,线性拟合所使用的基函数也未必一定是多项式, 根据实际问题可以设置为三角函数、指数函数、正态分布的概率密度函数,以及混合定义的基底函数。

  • 版权声明: 本博客所有文章除特别声明外,著作权归作者所有。转载请注明出处!

扫一扫,分享到微信

微信分享二维码
  • Copyrights © 2020-2021 chenk
  • 由 帅气的CK本尊 强力驱动
  • 访问人数: | 浏览次数:

请我喝杯咖啡吧~

支付宝
微信