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

您可以使用

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函数来获取从每个bootstrap样本计算的统计量,我们使用它来显示转换矩阵中每个条目的直方图：

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