matlab - マルコフ遷移行列の信頼区間の推定

Question

ACGTE という文字を含むさまざまな長さの一連のn=400シーケンスがあります。たとえば、A の後にCがある確率は次のとおりです。

そして、これは一連の経験的シーケンスから計算できるため、

仮定:

次に、遷移行列を取得します。

しかし、ファットの信頼区間を計算することに興味があります。どうすればそれを実行できるかについて何か考えがありますか?

Answer 1

ブートストラップを使用して信頼区間を推定できます。MATLAB はbootci、Statistics ツールボックスに機能を提供します。次に例を示します。

%# 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

および信頼区間の下限と上限を含む他の 2 つの同様の行列:

>> 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

Answer 2

統計的に正しいかどうかはわかりませんが、目安となる上限と下限を取得する簡単な方法は次のとおりです。

サンプルを n 等分 (1:40、41:80、...、361:400 など) に分割し、これらの各サブサンプルの確率行列を計算します。

サブサンプル間の確率の分布を見ることで、分散がどのようなものかをよく理解できるはずです。

この方法の欠点は、希望する特定の確率で間隔を実際に計算できない場合があることです。利点は、シリーズがどのように動作するかについて良い感触を与えること、および他の方法 (ブートストラップなど) が基づいているという仮定のために他の方法で失われる可能性のある情報を取得できることです。