了解 Matlab 中的伪随机数如何暗示统计独立性

Question

考虑以下 Matlab 代码，其中我使用伪随机数生成器生成一些数据。 我希望您能帮助我从统计的角度理解这些数字的“随机性”，我在下面解释。

我先设置了一些参数

%%%%%%%%Parameters
clear
rng default
Xsup=-1:6; 
Zsup=1:10; 
n_m=200; 
n_w=200; 
R=n_m;

然后我生成数据

%%%%%%%%Creation of data [XZ,etapair,zetapair,etasingle,zetasingle]

%Vector X of dimension n_mx1
idX=randi(size(Xsup,2),n_m,1); %n_mx1
X=Xsup(idX).'; %n_mx1

%Vector Z of dimension n_wx1
idZ=randi(size(Zsup,2),n_w,1); 
Z=Zsup(idZ).'; %n_wx1

%Combine X and Z in a matrix XZ of dimension (n_m*n_w)x2 
which lists all possible combinations of values in X and Z
[cX, cZ] = ndgrid(X,Z);
XZ = [cX(:), cZ(:)]; %(n_m*n_w)x2

%Vector etapair of dimension (n_m*n_w)x1
etapair=randn(n_m*n_w,1); %(n_m*n_w)x1

%Vector zetapair of dimension (n_m*n_w)x1
zetapair=randn(n_m*n_w,1); %(n_m*n_w)x1

%Vector etasingle of dimension (n_m*n_w)x1
etasingle=max(randn(n_m,R),[],2); %n_mx1 
etasingle=repmat(etasingle, n_w,1); %(n_m*n_w)x1

%Vector zetasingle of dimension (n_m*n_w)x1
zetasingle=max(randn(n_w,R),[],2); %n_wx1
zetasingle=kron(zetasingle, ones(n_m,1)); %(n_m*n_w)x1

现在让我将这些抽签转化为统计术语：

对于t=1,...,n_w*n_m，X(t)可以认为是随机变量X_t的实现

对于t=1,...,n_w*n_m，Z(t)可以认为是随机变量Z_t的实现

对于t=1,...,n_w*n_m，etapair(t)可以认为是随机变量E_t的实现

对于t=1,...,n_w*n_m，zetapair(t)可以认为是随机变量Q_t的实现

对于t=1,...,n_w*n_m，etasingle(t)可以认为是随机变量Y_t的实现

对于t=1,...,n_w*n_m，zetasingle(t)可以认为是随机变量S_t的实现

我的信念是，Matlab 中的伪随机数生成器允许声称 (X_1,X_2,..., Z_1,Z_2,...,E_1,E_2,..., Q_1,Q_2...,Y_1,Y_2,...,S_1,S_2,...) 相互独立 如解释here

作为对这个假设声明的检查，我定义了 W_t:=-E_t-Q_t+Y_t+S_t 并根据经验计算 Pr(W_t<=1|X_t=5, Z_t=1)

如果相互独立成立，那么Pr(W_t<=1|X_t=5, Z_t=1)=Pr(W_t<=1) 及其下面的经验对应物，命名为option1 和option2，应该几乎相同。

%option 1
num1=zeros(n_m*n_w,1);
for h=1:n_m*n_w
    if -etapair(h)-zetapair(h)+etasingle(h)+zetasingle(h)<=1 && XZ(h,1)==5 && XZ(h,2)==1
        num1(h)=1;
    end
end
den1=zeros(n_m*n_w,1);
for h=1:n_m*n_w
    if  XZ(h,1)==5 && XZ(h,2)==1
        den1(h)=1;
    end
end
option1=sum(num1)/sum(den1);

%option 2
num2=zeros(n_m*n_w,1);
for h=1:n_m*n_w
    if -etapair(h)-zetapair(h)+etasingle(h)+zetasingle(h)<=1 
        num2(h)=1;
    end
end
option2=sum(num2)/(n_m*n_w);

问题：option1 (=0.0021) 和 option2 (=0.0012) 之间的区别是指“几乎”还是我做错了什么？

Answer 1

由于观察随机事件的本质，你不能保证给定经验试验的理论上准确的结果。

您在脚本的开头设置了rng default，这意味着您将始终得到相同的结果（option1 = 0.0021、option2 = 0.0012）。

多次运行脚本并对结果取平均值，我们应该接近理论准确性。

kk = 10000;
option1 = zeros(kk, 1);
option2 = zeros(kk, 1);
for ii = 1:kk
    % No need to use 'clear' here. If you were concerned 
    % for some reason, you could use 'clearvars -except kk option1 option2 ii'
    % do not use 'rng default'. Use 'rng shuffle' if anything, but not necessary
    Xsup = -1:6;
    % ... all your other code
    % replace 'option1=...' with 'option1(ii)=...'
    % replace 'option2=...' with 'option2(ii)=...'  
end
fprintf('Results:\nMean option1 = %f\nMean option2 = %f\n', mean(option1), mean(option2));

结果：

>> Mean option1 = 0.001461
>> Mean option2 = 0.001458

我们可以看到这些在一定程度上是相同的，如果我们运行 X 次试验（对于足够大的 X），准确度可以任意高。这与自变量的预期一致。

---

注意，如果你有并行计算工具箱，这个for 循环可以很容易地换成parfor，你可以更快地运行试验。