matlab - 估计马尔可夫转移矩阵的置信区间

Question

我有一系列不同长度的 n=400 序列，其中包含字母 ACGTE。例如，在 A 之后有 C 的概率是:

并且可以从经验序列的集合中计算出，因此

假设:

然后我得到一个转换矩阵:

但我对计算 Phat 的置信区间很感兴趣，有什么想法可以解决吗？

Answer 1

你可以使用 bootstrapping估价confidence intervals . MATLAB 提供 bootci统计工具箱中的函数。这是一个例子:

%# generate a random cell array of 400 sequences of varying length
%# each containing indices from 1 to 5 corresponding to ACGTE
sequences = arrayfun(@(~) randi([1 5], [1 randi([500 1000])]), 1:400, ...
    'UniformOutput',false)';

%# compute transition matrix from all sequences
trans = countFcn(sequences);

%# number of bootstrap samples to draw
Nboot = 1000;

%# estimate 95% confidence interval using bootstrapping
ci = bootci(Nboot, {@countFcn, sequences}, 'alpha',0.05);
ci = permute(ci, [2 3 1]);

我们得到:

>> trans         %# 5x5 transition matrix: P_hat
trans =
      0.19747       0.2019      0.19849       0.2049      0.19724
      0.20068      0.19959      0.19811      0.20233      0.19928
      0.19841      0.19798       0.2021       0.2012      0.20031
      0.20077      0.19926      0.20084      0.19988      0.19926
      0.19895      0.19915      0.19963      0.20139      0.20088

和其他两个包含置信区间下限和上限的类似矩阵:

>> ci(:,:,1)     %# CI lower bound
>> ci(:,:,2)     %# CI upper bound

我正在使用以下函数从一组序列中计算转换矩阵:

function trans = countFcn(seqs)
    %# accumulate transition matrix from all sequences
    trans = zeros(5,5);
    for i=1:numel(seqs)
        trans = trans + sparse(seqs{i}(1:end-1), seqs{i}(2:end), 1, 5,5);
    end

    %# normalize into proper probabilities
    trans = bsxfun(@rdivide, trans, sum(trans,2));
end

作为奖励，我们可以使用 bootstrp 函数来获取从每个引导样本计算的统计数据，我们用它来显示转换矩阵中每个条目的直方图:

%# compute multiple transition matrices using bootstrapping
stat = bootstrp(Nboot, @countFcn, sequences);

%# display histogram for each entry in the transition matrix
sub = reshape(1:5*5,5,5);
figure
for i=1:size(stat,2)
    subplot(5,5,sub(i))
    hist(stat(:,i))
end