Implementating the use of strfind with a for loop

I am attempting to write a transition matrix and to do so I will be comparing the current state to the next state. It seems like using strfind and a couple for loops would be the best way to do this, but I don't know how to. Any suggestions on how to begin?

Thanks

[James Tursa]

Please give a short example of the input variables and a result you want, with a description of what the algorithm is supposed to do.

Answer 1

[arich82]

[Edited to correct sum(p, 1) to sum(p, 2).]

I think you submitted a similar question yesterday, but it seems to have been deleted. (I saw it, but the forum gave me an error when I tried to log on).

Warning: I know virtually nothing of Markov chains and transition matrices (so naturally, I'm going to try to address the question)...

You'd posted the example data:

d = [3 3 4 4 4 5 6 5 5 4];

and I think you want a matrix p of the probability that a given entry is followed by a given value.

[Note: if you want to use strfind like you did yesterday, I think you need to look at the help doc: it returns a list of of the first index of the start of the substring, so you'd want to use numel on the result to find how many substrings there were, not sum (which would sum the values of the indices, not their number).]

It makes more sense to me to use the additional output from the unique function that you used previously, then follow the method described here.

d = [3 3 4 4 4 5 6 5 5 4];
[c, ~, ic] = unique(d); 
% c is a vector of the unique values in d, while 
% ic is the same length as d, and is the index of c 
% that gives the original value of d, i.e. d = c(ic);
% this allows us to index directly into p, even if our 
% data is non-consecutive, or even non-integral
% (essentially, ic converts the values in d to arbitrary integer labels)

N = numel(c); % number of unique entries in d
p = zeros(N, N); % preallocate transition matrix

% cycle through the data once, looking at the current and next value, 
% i.e. ic(k) and ic(k + 1), using p to tally the number of times 
% a given combination was encountered
for k = 1:(numel(d) - 1)
  p(ic(k), ic(k + 1)) = p(ic(k), ic(k + 1)) + 1;
end

% normalize the count stored in p by the number of times a given
% starting value was encountered, i.e. the sum of each row
p = bsxfun(@rdivide, p, sum(p, 2));

% Note: the preceding bsxfun command is equivalent to the loop:
%   for k = 1:N
%     p(k, :) = p(k, :) / sum(p(k, :));
%   end
%

The result is

p =

      0.5000    0.5000         0         0
           0    0.6667    0.3333         0
           0    0.3333    0.3333    0.3333
           0         0    1.0000         0

I interpret the first row to mean: given a 3, there is a 50% chance the next value will be another 3, a 50% chance the next value will be a 4, and 0% chance of anything else. I interpret the second row to mean: given a 4, there is a 67% chance the next value will be another 4, a 33% chance the next value will be a 5, and a 0% chance of anything else. An so on... (Is this what you want?)

If you really insist on using strfind (and assuming the above procedure computed the correct result), I think the following should be equivalent (but far less efficient since strfind will cycle through the data N^2 times, instead of only once), and is closer to what you started with yesterday:

d = [3 3 4 4 4 5 6 5 5 4];
u1 = unique(d);
N = numel(u1);
x = zeros(N, 1);
y = zeros(N, N);
for i = 1:N
    for j = 1:N
        y(i, j) = numel(strfind(d, [u1(i), u1(j)]));
    end
    x(i) = sum(d(1:end - 1) == u1(i));
end
p = bsxfun(@rdivide, y, x);

% Note: Matlab uses Fortran (column-based) storage for arrays,
% so nested loops are usually ordered j, then i;
% also, using i and j for indices is sometimes discouraged, since
% they are used for complex number by default, though it's ok if you
% **always** use 1i to express an imaginary number insted a 1*i;
% that said, I've re-used your code from yesterday as much as possible.
%
% Finally, this is a terribly inefficient way of getting the count in x;
% note that you also must exclude the last value in d for the count,
% since you don't have stats on what follows it

Please accept this answer if it helps, or let me know in the comments if it completely misses the mark...

--Andy

[arich82]

After thinking about it, I realized this is essentially just a 2-D binning operation, where the pairs are just [ic(k), ic(k+1)].

The following code is a little opaque, but should be very fast, and doesn't require any toolboxes:

d = [3 3 4 4 4 5 6 5 5 4];
[c, ~, ic] = unique(d);
N = numel(c);
p = accumarray([ic(1:end - 1), ic(2:end)], 1, [N, N]);
p = bsxfun(@rdivide, p, sum(p, 2));

Not sure if it's helpful, but I'll leave it here for posterity...

[Ellie]

Is there any way that I could modify this code so that it is [[ic(k), ic(k+1), ic(k+2). ic(k+3)]?

Answer 2

[Steven Lord]

You may find the Hidden Markov Models functionality in Statistics and Machine Learning Toolbox (previously called Statistics Toolbox) useful.