18、MATLAB中的随机数、分布、映射与噪声模拟

MATLAB中的随机数、分布、映射与噪声模拟

1. 随机数与分布

1.1 MATLAB中的协方差和相关系数命令

MATLAB提供了 cov 和 corrcoef 命令。 cov(x,y) 用于计算两个向量之间的协方差，并返回一个包含方差的矩阵：
[
cov(x,y) =
\begin{bmatrix}
\sigma_X & \sigma_{XY} \
\sigma_{YX} & \sigma_Y
\end{bmatrix}
]
corrcoef(x,y) 则用于计算相关系数矩阵：
[
corrcoef(x,y) =
\begin{bmatrix}
1 & r_{XY} \
r_{YX} & 1
\end{bmatrix}
]
对角线上的元素表明随机变量与其自身是相关的。

1.2 随机数生成

1.2.1 rand 命令

rand 命令可生成0到1之间的随机数。以下是一个简单的掷骰子程序示例：

% roll.m
function [number] = roll()
number = ceil(rand*6+eps);

调用该函数的代码如下：

% mroll.m
for ii = 1:6000
    rr(ii) = roll;
end
for j = 1:6
    ii(j) = length(find(rr==j));
end
disp([' Mean ' num2str(mean(rr))])

这里加入 eps 是为了防止随机数为0的情况。样本均值的计算公式为：
[
\bar{x} = \frac{1}{N} \sum_{i=1}^{N} x_i
]
理论上掷骰子的均值为7/2。样本方差的计算公式为：
[
var = \frac{1}{N - 1} \sum_{i=1}^{N} (x_i - \bar{x})^2
]
对于掷骰子的情况，方差约为35/12 ≈ 2.9167。

如果想模拟两个骰子的情况，只需将 mroll.m 中的 rr(ii) = roll; 修改为 rr(ii) = roll+roll; 。

1.2.2 抛硬币模拟

以下是一个抛硬币的模拟函数：

function [out] = toss()
A = ['heads';'tails'];
i = ceil(rand*2+eps);
out = A(i,:);

调用 toss 函数即可生成 heads 或 tails 。

1.3 正态分布

正态分布有两个参数：均值和方差。MATLAB中的 randn 命令可生成均值为0、方差为1的随机数。示例代码如下：

x = randn(1,500);
hist(x,12);

可以使用 mean 和 var 命令来分析生成的随机数与理论值的接近程度：

>> mean(x)
ans =
    0.0702
>> var(x)
ans =
    0.9873

1.3.1 示例：分析工资分布

考虑工资分布，均值为£20000，标准差为£3000。要计算200人中工资超过£22000的人数，可使用以下代码：

no_experiments = 200;
for ii = 1:no_experiments
    x = randn(1,200);
    salaries = x*3000+20000;
    number(ii) = length(find(salaries>22000));
end
mean(number)

运行结果约为50.94。

1.4 计算概率

随机变量X小于值z的概率计算公式为：
[
P(X \leq z) = \frac{1}{\sigma \sqrt{2\pi}} \int_{-\infty}^{z} e^{-\frac{1}{2}(\frac{q - \mu}{\sigma})^2} dq
]
以下是使用MATLAB计算该概率的代码：

mu = 0.2;
sig = 1.5;
min_inf = -6*sig;
z = input('Enter z: ');
clear x n
x = linspace(min_inf,z,1000);
n = length(x);
h = x(2)-x(1);
f = exp(-((x-mu)/sig).^2/2);
val = (sum(f)-f(1)/2-f(n)/2)*h;
prob = val/sqrt(2*pi)/sig

该代码使用了梯形法则，并将积分下限设置为 -6σ。

1.5 排列

MATLAB提供了 randperm 命令用于生成随机排列。例如，要对五个名字进行随机排序：

a = ['Avril'; 'Beryl'; 'Carol'; 'Daisy'; 'Emily'];
r = randperm(5);
a(r,:)

2. 映射与白噪声

2.1 一维映射

一维映射可以表示为 ( x_{n+1} = f(x_n) )。固定点满足 ( x_n = f(x_n) )。例如，帐篷映射的定义如下：
[
f(x) =
\begin{cases}
ax, & x \leq \frac{1}{2} \
a(1 - x), & x > \frac{1}{2}
\end{cases}
]
可改写为 ( f(x) = a(\frac{1}{2} - |x - \frac{1}{2}|) )。以下是实现该映射的MATLAB代码：

function [out] = tent(a,x)
out = a*(0.5-abs(0.5-x));

绘制 ( a = 4 ) 时的映射结果：

x = linspace(-1,2);
plot(x,tent(4,x))
h = gca;
set(h,'YGrid','On')

2.2 二维映射

考虑一个简单的二维映射：
[
\begin{cases}
\dot{x} = y \
\dot{y} = -x
\end{cases}
]
使用Crank - Nicolson方法离散化该系统：
[
\begin{bmatrix}
1 & -\frac{\Delta t}{2} \
\frac{\Delta t}{2} & 1
\end{bmatrix}
\begin{bmatrix}
x_{n+1} \
y_{n+1}
\end{bmatrix}
=
\begin{bmatrix}
x_n + \frac{\Delta t}{2} y_n \
y_n - \frac{\Delta t}{2} x_n
\end{bmatrix}
]
以下是实现该映射的MATLAB代码：

% maps.m
function [f,g]=maps(x,y)
global delt
A = [ 1 -delt/2; delt/2 1];
rhs = [x+delt/2*y; -delt/2*x+y];
v = A\rhs;
f = v(1);
g = v(2);

% crmap.m
global delt
delt = 0.05;
x = 0.5; y = 0.4;
noise = 0.0;
nop = 500;
po(1) = x+i*y;
for it = 2:nop
    [xn,yn] = maps(x,y);
    xn = xn + noise*sqrt(-2*log(rand(1)))*cos(2*pi*rand(1));
    yn = yn + noise*sqrt(-2*log(rand(1)))*sin(2*pi*rand(1));
    po(it) = xn+i*yn;
    x = xn; y = yn;
end
plot(po,'.')
axis equal

该系统的解应该是一个圆。

2.3 离散系统建模

考虑一个捕食者 - 猎物模型：
[
\begin{cases}
N_{t+1} = rN_t e^{-aP_t} \
P_{t+1} = N_t (1 - e^{-aP_t})
\end{cases}
]
以下是求解该模型的MATLAB代码：

r = 2;
a = 1;
nt = 20;
t = 1:1:nt;
N = zeros(size(t));
P = zeros(size(t));
N(1) = r*log(r)/(a*(r-1))+0.01;
P(1) = log(r)/a;
for j = 1:(nt-1)
    N(j+1) = r*N(j)*exp(-a*P(j));
    P(j+1) = N(j)*(1-exp(-a*P(j)));
end
subplot(2,1,1)
plot(t,N,'o','MarkerSize',12)
title('N values')
subplot(2,1,2)
plot(t,P,'r+','MarkerSize',12)
title('P values')

结果显示，稍微改变平衡振幅的值，种群会开始振荡并最终偏离初始值，说明该平衡是不稳定的。

2.4 周期性与混沌

考虑一个种群模型：
[
N_{t+1} = N_t + rN_t(1 - N_t)
]
该差分方程的解取决于参数r的值。以下是求解该模型的MATLAB代码：

str = 'Please enter the ';
r = input([str 'rate constant r :']);
N0 = input([str 'initial population N(0) :']);
n = input([str 'number of years n :']);
t = 1:1:n;
N = zeros(size(t));
N(1) = N0;
for j = 1:(n-1)
    N(j+1) = N(j) + r*N(j)*(1-N(j));
end
plot(t-1,N,'*','MarkerSize',10)

当r取不同值时，系统会表现出不同的行为：
- 当r较小时，种群趋于稳定值。
- 当r增大时，会出现周期振荡。
- 当r进一步增大时，会出现混沌现象。

2.5 随机运动

考虑一个粒子在盒子中的运动，西部墙壁对粒子有排斥力，东部墙壁位于x = 5。粒子的位置更新公式为：
[
x_{t+1} = x_t + f(x_t) + Noise
]
其中 ( f(x) = \frac{k}{x^2} )，噪声为 ( \epsilon ) 乘以 -1 到 1 之间的随机数。以下是模拟该运动的MATLAB代码：

x0 = 2.5;
noise = 1.;
east_wall = 5;
k = 0.5;
no_experiments = 200;
for j = 1:no_experiments
    x = x0;
    steps = 0;
    while x<east_wall
        steps=steps+1;
        x = x + k./x^2 + noise*(rand(1)-0.5)*2;
    end
    st_store(j) = steps;
end

可以使用 mean 、 std 和 hist 命令对实验结果进行统计分析。

2.6 数值求解周期点

2.6.1 迭代映射函数

为了找到特定周期的点，可以编写一个函数来返回点经过n次迭代后的图像。以下是实现该功能的MATLAB代码：

% maplog
function [new] = maplog(old,r,n)
for j=1:n
    new = old+r*old*(1-old);
    old = new;
end

2.6.2 牛顿 - 拉夫森方法求解

使用牛顿 - 拉夫森方法来数值求解周期点。以下是主程序代码：

r = input('Please enter the value of r :');
starting_guess = input('Please enter starting guess ');
ms = 'Please enter the period of the orbit sought ';
period = input(ms);
maxits = 200;
tolerance = 1e-10;
delta = 1e-4;
x = starting_guess;
for its = 1:maxits
    err = maplog(x,r,period)-x;
    if abs(err)<tolerance
        break
    end
    err_dx = maplog(x+delta,r,period)-(x+delta);
    x = x - delta*err/(err_dx-err);
end
if its==maxits
    disp(['No period ' int2str(period) ' point found ']);
else
    disp('Orbit');
    for j = 1:period
        disp(x);
        x = maplog(x,r,1);
    end
end

例如，当r = 2.1，周期为2时，输入合适的初始猜测值，程序会输出近似的周期点，约为1.1286和0.8237。

2.6.3 绘制固定点图

可以绘制该映射的固定点图，展示不同r值下系统的行为。以下是绘制固定点图的MATLAB代码：

clf
for r = 0:0.01:4;
    x = rand(1);
    % Remove transients
    for j = 1:100
        x = x+r*x*(1-x);
    end
    xout = [];
    for j = 1:400
        x = x+r*x*(1-x);
        xout = [xout x];
    end
    plot(r*ones(size(xout)),xout,'.','MarkerSize',3);
    hold on
end
axis([0 4 0 1])

运行该代码后，可以得到一个展示不同r值下系统固定点的图形。当放大图形的特定区域时，可以更清晰地观察到系统的分岔和混沌现象。例如，在r值逐渐增大的过程中，系统会经历周期加倍分岔，最终进入混沌状态。

3. 总结

3.1 主要内容回顾

本文主要介绍了MATLAB在随机数生成、分布模拟、映射建模以及噪声分析等方面的应用。具体内容如下：
| 主题 | 详细内容 |
| ---- | ---- |
| 随机数与分布 | 介绍了 rand 、 randn 等命令生成随机数，涵盖均匀分布、正态分布等。还包括计算协方差、相关系数、概率等方法，以及抛硬币、掷骰子等模拟示例。 |
| 映射与白噪声 | 探讨了一维映射（如帐篷映射）和二维映射的建模与分析，使用Crank - Nicolson方法离散化系统。还介绍了离散系统建模（捕食者 - 猎物模型）、周期性与混沌（种群模型）以及随机运动（粒子在盒子中的运动）的模拟。 |
| 数值求解 | 介绍了使用牛顿 - 拉夫森方法数值求解映射的周期点，并绘制固定点图展示系统的行为。 |

3.2 操作流程总结

以下是本文涉及的主要操作流程：
1. 随机数生成与分布模拟
- 使用 rand 生成0到1之间的随机数，用于模拟掷骰子、抛硬币等。
- 使用 randn 生成均值为0、方差为1的正态分布随机数。
- 使用 cov 和 corrcoef 计算协方差和相关系数。
- 使用积分方法计算随机变量的概率。
2. 映射建模与分析
- 定义一维和二维映射函数，使用迭代方法计算映射的后续值。
- 使用Crank - Nicolson方法离散化二维系统。
- 分析映射的固定点、周期点和混沌行为。
3. 数值求解与绘图
- 使用牛顿 - 拉夫森方法数值求解映射的周期点。
- 绘制映射的固定点图，观察系统的分岔和混沌现象。

3.3 整体流程mermaid流程图

graph LR
    classDef startend fill:#F5EBFF,stroke:#BE8FED,stroke-width:2px;
    classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;
    classDef decision fill:#FFF6CC,stroke:#FFBC52,stroke-width:2px;

    A([开始]):::startend --> B(随机数与分布模拟):::process
    B --> B1(生成随机数):::process
    B --> B2(计算协方差和相关系数):::process
    B --> B3(计算概率):::process
    B --> B4(模拟实验):::process
    A --> C(映射与白噪声建模):::process
    C --> C1(定义映射函数):::process
    C --> C2(离散化系统):::process
    C --> C3(分析固定点和周期点):::process
    C --> C4(观察混沌行为):::process
    A --> D(数值求解与绘图):::process
    D --> D1(数值求解周期点):::process
    D --> D2(绘制固定点图):::process
    B --> E{是否完成所有实验?}:::decision
    C --> E
    D --> E
    E -->|否| F(继续实验或分析):::process
    F --> B
    F --> C
    F --> D
    E -->|是| G([结束]):::startend

通过以上内容，我们可以看到MATLAB在处理随机现象和动态系统建模方面具有强大的功能。合理运用这些功能，可以帮助我们更好地理解和分析各种复杂的系统。希望本文的介绍对读者在相关领域的研究和实践有所帮助。