如何计算矩阵中一个值的连续出现长度（正常运行时间）？

Question

我有这样的数据:

 1     0     1
 1     1     1
 0     1     1
 1     1     1
 1     1     1
 1     1     1
 1     1     0
 1     1     1
 1     1     1
 1     1     1
 1     1     1
 1     1     1
 1     1     1
 1     1     1
 0     0     1
 1     1     1
 1     1     1
 1     1     1

每一列代表一个设备，每一行代表一个时间段。每个数据点表示该时间段内设备是否处于活动状态。我正在尝试计算每个设备连续处于活动状态的时间长度，或称为“连续时间”，即在每一列中连续的 1 的长度。在这种情况下，第一列应该是2 11 3，依此类推。

对于单个设备（单列数据），这很容易实现：

rng(1)

%% Parameters
lambda = 0.05;      % Pr(failure)
N = 1;              % number of devices
T = 18;             % number of time periods in sample

%% Generate example data
device_status = [rand(T, N) >= lambda ; false(1, N)];

%% Calculate spell lengths, i.e. duration of uptime for each device
cumul_status = cumsum(device_status);

% The 'cumul_status > 0' condition excludes the case where the vector begins with one
% or more zeros
cumul_uptimes = cumul_status(device_status == 0 & cumul_status > 0);
uptimes = cumul_uptimes - [0 ; cumul_uptimes(1:end-1)];

所以我可以简单地迭代遍历列，并逐一进行操作，使用parfor（例如）并行运行。是否有一种方法可以使用向量化矩阵操作同时在所有列上执行此操作？

编辑：应该补充说明的是，每个设备可能具有不同数量的正常运行时间。

- Michael A

---

你可以尝试使用类似以下的代码：reshape(max(cumsum(cumprod(hankel([zeros(1,size(A,2));A])'))),size(A,1)+1,size(A,2)); 其中 A 是矩阵。然后你就可以找到局部最大值，但有时使用循环更好更快。 - obchardon

如果您没有太多的设备，那么我认为使用循环的这种方法是完全可以的。不要过于复杂化事情。Donald Knuth表达得最好：“过早优化是万恶之源”。除非您可以清楚地看到逐个循环遍历每个设备与试图以矢量化方式处理所有设备之间存在巨大的运行时差异，否则这样做并不值得。 - rayryeng

另外需要注意的是，每列的正常运行时间可能不同。具体来说，在您的示例中，第一列有三个正常运行时间实例，而最后两列只有两个。这种“不平衡”已经影响了您进行向量化的能力。 - rayryeng

@rayryeng 非常好的观点。我一开始就知道这一点，但是我忘记在问题中清楚地说明了。我添加了一个编辑。此外，你说循环可能是最好的选择，你是对的。我主要是为了教育而提出这个问题，因为每当我在SO上提问时，我总是学到很多东西。 - Michael A

Answer 1

这里有一种方法，不确定是否算作向量化。

将数据矩阵表示为x。

[ii, jj] = find([true(1,size(x,2)); ~x; true(1,size(x,2))]);
result = accumarray(jj, ii, [], @(x){nonzeros(diff(x)-1)});

该函数生成一个单元数组，其中每个单元格对应一列。在你的例子中，

result{1} =
     2
    11
     3
result{2} =
    13
     3
result{3} =
     6
    11

如何运作 这个想法是找到x中零的行和列索引（也就是在~x中为true的值），然后使用列索引作为分组变量（accumarray的第一个参数）。

在每个组内，我们使用匿名函数@(x){nonzeros(diff(x)-1)}来计算零的行位置之间的差异。我们可以直接应用diff，因为从find得到的列索引已经排序好了，这要归功于Matlab的列优先顺序。我们减去1，因为x中的零不算作正常运行时间；删除长度等于0的正常运行时间（使用nonzeros），并将结果打包成一个单元格（{...}）。

在~x的前后添加一行true值，以确保我们检测到初始和最终的正常运行时间段。

- Luis Mendo