6

我有一个以下随机模型，描述了一个过程（Y）在空间和时间上的演化。Ds 和 Dt 是空间域（具有 x 和 y 轴的 2D）和时间域（具有 t 轴的 1D）。该模型通常称为混合效应模型或变异分量模型

在此处输入图像描述

我目前正在开发 Y 如下：

%# Time parameters
T=1:1:20; % input
nT=numel(T);

%# Grid and model parameters
nRow=100;
nCol=100;


[Grid.Nx,Grid.Ny,Grid.Nt] = meshgrid(1:1:nCol,1:1:nRow,T);

xPower=0.1;
tPower=1;
noisePower=1;
detConstant=1;

deterministic_mu = detConstant.*(((Grid.Nt).^tPower)./((Grid.Nx).^xPower));

beta_s = randn(nRow,nCol); % mean-zero random effect representing location specific variability common to all times

gammaTemp = randn(nT,1);

for t = 1:nT
    gamma_t(:,:,t) = repmat(gammaTemp(t),nRow,nCol); % mean-zero random effect representing time specific variability common to all locations
end

var=0.1;% noise has variance = 0.1
for t=1:nT
    kappa_st(:,:,t) = sqrt(var)*randn(nRow,nCol);
end

for t=1:nT
    Y(:,:,t) = deterministic_mu(:,:,t) + beta_s + gamma_t(:,:,t) + kappa_st(:,:,t);
end

我的问题是：

  1. 如何在 Y 的表达式中产生 delta 以及 kappa 和 delta 的差异？
  2. 通过使用 Matlab 的一些说明帮助解释我是否正确生成 Y？

如果您需要更多信息/解释，请告诉我。谢谢。

4

2 回答 2

4

首先，我重写了您的代码以使其更有效率。我看到您为 生成线性间隔网格x，y并t为此网格中的所有点执行计算。这种方法对可达到的最大网格分辨率有严格的限制，因为如果分辨率提高，3D 网格（以及用它定义的所有变量）会消耗大量内存。如果您正在实现的模型的复杂性和大小会增加（它经常如此），我建议您将这一切都投入到一个接受矩阵/向量输入的函数中s，t这在这方面会更加灵活——以这种方式处理原本无法放入内存的数据“块”会容易得多。

然后，我生成了这个delta_st术语，rand而不是randn因为噪声应该是“白色”。现在我对最后一个非常不确定，而且我没有时间阅读你链接到的论文——你能告诉我在哪些页面上我可以找到相关的部分delta_st吗？

现在，代码：

%# Time parameters
T  = 1:1:20; % input
nT = numel(T);

%# Grid and model parameters
nRow = 100;
nCol = 100;

% noise has variance = 0.1
var = 0.1;

xPower = 0.1;
tPower = 1;
noisePower  = 1;
detConstant = 1;



[Grid.Nx,Grid.Ny,Grid.Nt] = meshgrid(1:nCol,1:nRow,T);

% deterministic mean
deterministic_mu = detConstant .* Grid.Nt.^tPower ./ Grid.Nx.^xPower;

% mean-zero random effect representing location specific 
% variability common to all times
beta_s = repmat(randn(nRow,nCol), [1 1 nT]);

% mean-zero random effect representing time specific 
% variability common to all locations
gamma_t = bsxfun(@times, ones(nRow,nCol,nT), randn(1, 1, nT));


% mean zero random effect capturing the spatio-temporal 
% interaction not found in the larger-scale deterministic mu 
kappa_st = sqrt(var)*randn(nRow,nCol,nT);


% mean zero random effect representing the micro-scale
% spatio-temporal variability that is modelled by white
% noise (i.i.d. at different time steps) in Ds·Dt
delta_st = noisePower * (rand(nRow,nCol,nT)-0.5);


% Final result: 
Y = deterministic_mu + beta_s + gamma_t + kappa_st + delta_st;
于 2012-10-02T13:02:55.590 回答
1

您的实现对 beta、gamma 和 kappa 进行采样，就好像它们是白色的一样（例如，它们在每个 (x,y,t) 处的值是独立的）。这些术语的描述表明并非如此。看起来 delta 应该捕获白噪声，而其他项捕获它们各自域上的相关性。例如，gamma(t_1) 和 gamma(t_1+1) 之间存在非零相关性。

如果您希望将 gamma 建模为具有方差 var_g 和 gamma(t) 和 gamma(t+1) 之间的相关 cor_g 的平稳高斯马尔可夫过程，您可以使用类似

gamma_t = nan( nT, 1 );
gamma_t(1) = sqrt(var_g)*randn();
K_g = cor_g/var_g;
K_w = sqrt( (1-K_g^2)*var_g );
for t = 2:nT,
   gamma_t(t) = K_g*gamma_t(t-1) + K_w*randn();
end
gamma_t = reshape( gamma_t, [ 1 1 nT ] );

我在上述代码中用于增益 K_g 和 K_w 的公式（以及 gamma_t(1) 的初始化）产生了所需的平稳方差 \sigma^2_0 和一步协方差 \sigma^2_1：

增益公式，并证明方差和协方差符合要求

请注意，上面的实现假设稍后您将使用 bsxfun 为您执行“repmat”对术语求和：

Y = bsxfun( @plus, deterministic_mu + kappa_st + delta_st, beta_s );
Y = bsxfun( @plus, Y, gamma_t );

请注意，我没有测试过上面的代码，因此您应该通过采样确认它确实会产生相邻样本之间指定方差和协方差的零噪声过程。对 beta 进行采样，同样的过程可以扩展到二维，但原理基本相同。我怀疑 kappa 应该被类似地建模为马尔可夫高斯过程，但在所有三个维度上并具有较低的方差来表示未在 mu、beta 和 gamma 中捕获的高阶效应。

Delta 应该是零均值静止白噪声。假设它是高斯方差，noisePower 将使用它对其进行采样

delta_st = sqrt(noisePower)*randn( [ nRows nCols nT ] );
于 2012-10-06T08:50:52.837 回答