Matlab——带边界的随机游走，矢量化

Question

假设我有一个跳跃大小的向量 J 和一个初始起点 X_0。我也有边界 0，B（假设 0

在下面的代码中，我在许多示例中执行此操作。要“修复”超出边界的那些，我必须遍历样本以检查...（不要认为有矢量化的“查找”）

% X_init is a row vector describing initial resource values to use for
% each sample

% J is matrix where each col is a sequence of Jumps (columns = sample #)
% In this code the jumps are subtracted, but same thing

X_intvl = repmat(X_init,NumJumps,1) - cumsum(J);
X = [X_init; X_intvl];
for sample = 1:NumSamples
    k = find(or(X_intvl(:,sample) > B, X_intvl(:,sample) < 0),1);
    while(~isempty(k))
        change = X_intvl(k-1,sample) - X_intvl(k,sample);
        X_intvl(k:end,sample) = X_intvl(k:end,sample)+change;
        k = find(or(X_intvl(:,sample) > B, X_intvl(:,sample) < 0),1);
    end
end

Answer 1

有趣的问题 (+1)。

前段时间我也遇到过类似的问题，虽然由于我的下限和上限取决于 t，所以稍微复杂一些。我从来没有制定出完全矢量化的解决方案。最后，我发现最快的解决方案是单个循环，它在每一步都包含约束。根据您的情况调整代码会产生以下结果：

%# Set the parameters
LB = 0; %# Lower bound
UB = 5; %# Upper bound
T = 100; %# Number of observations
N = 3; %# Number of samples
X0 = (1/2) * (LB + UB); %# Arbitrary start point halfway between LB and UB

%# Generate the jumps
Jump = randn(N, T-1);

%# Build the constrained random walk
X = X0 * ones(N, T);
for t = 2:T
    X(:, t) = max(min(X(:, t-1) + Jump(:, t-1), UB), 0);
end
X = X';

我很想知道这种方法是否比您目前正在做的更快。我怀疑这将适用于约束在不止一个或两个地方具有约束力的情况。我无法自己测试它，因为您提供的代码不是“工作”示例，即我不能只是将其复制并粘贴到 Matlab 中并运行它，因为它取决于示例（或模拟）值的几个变量不提供。我尝试自己调整它，但无法使其正常工作？

更新：我刚刚切换了代码，以便观察在列上被索引，样本在行上被索引，然后我在最后一步中转置X。这将使例程更高效，因为 Matlab 为数字数组按列分配内存 - 因此在沿数组的列（而不是跨行）执行操作时更快。请注意，您只会注意到大型 N 的加速。

最后的思考：如今，JIT 加速器非常擅长在 Matlab 中提高单循环的效率（双循环仍然很慢）。因此，我个人认为，每次您尝试在 Matlab 中获得一个完全矢量化的解决方案时，即没有循环，您应该权衡寻找一个聪明的解决方案所涉及的努力是否值得稍微提高效率通过使用单个循环的更容易获得的方法。重要的是要记住，当 T 和 N 很小时，完全矢量化的解决方案有时会慢比涉及单循环的解决方案！

Answer 2

我想提出另一种矢量化解决方案。

所以，首先我们应该设置参数并生成随机Jumpls。我使用了与 Colin T Bowers 相同的一组参数：

% Set the parameters
LB = 0; % Lower bound
UB = 20; % Upper bound
T = 1000; % Number of observations
N = 3; % Number of samples
X0 = (1/2) * (UB + LB); % Arbitrary start point halfway between LB and UB

% Generate the jumps
Jump = randn(N, T-1);

但我更改了生成代码：

% Generate initial data without bounds
X = cumsum(Jump, 2);

% Apply bounds
Amplitude = UB - LB;
nsteps = ceil( max(abs(X(:))) / Amplitude - 0.5 );
for ii = 1:nsteps
    ind = abs(X) > (1/2) * Amplitude;
    X(ind) = Amplitude * sign(X(ind)) - X(ind);
end

% Shifting X
X = X0 + X;

因此，我使用带有智能后处理功能的 cumsum function 而不是 for 循环。

注意此解决方案的工作速度明显慢于 Colin T Bowers 的解决方案，适用于严格界限 (Amplitude < 5)，但适用于松散界限 (Amplitude > 20)更快。