CODEWITHSUNDEEP
matlabmarkov-chains
Ram s
Ram s

Reputation: 21

how to programme in Matlab to obtain transition probability matrix?

I have a sequence x= [12,14,6,15,15,15,15,6,8,8,18,18,14,14] so I want to make transition probability matrix. Transition probability matrix calculated by equation i.e. probability=(number of pairs x(t) followed by x(t+1))/(number of pairs x(t) followed by any state). Matrix should be like below

                  6      8    12     14     15      18

           6      0     1/2   0       0     1/2     0
           8      0     1/2   0       0      0     1/2
           12     0      0    0       1      0      0
           14    1/2     0    0      1/2     0      0
           15    1/4     0    0       0     3/4     0
           18     0      0    0       0     1/2    1/2

by following code I can do

m = max(x);
n = numel(x);
y = zeros(m,1);
p = zeros(m,m);
for k=1:n-1
    y(x(k)) = y(x(k)) + 1;
    p(x(k),x(k+1)) = p(x(k),x(k+1)) + 1;
end
p = bsxfun(@rdivide,p,y); p(isnan(p)) = 0;

but with this code matrix forms of order maximum state present in sequence i.e. matrix becomes of 18*18, and much more places zero occurs. I want matrix like above posted by me how to do it.

Upvotes: 0

Views: 1059

Answers (3)

Muhammad Zahid
Muhammad Zahid

Reputation: 11

x=[12,14,6,15,15,15,15,6,8,8,18,18,14,14]; %discretized driving cycle
n=numel(x);%total no of data points in driving cycle
j=0;
z=unique(x); % unique data points in the driving cycle and also in arranged form
m=numel(z); % total number of unique data points
N=zeros(m); % square matrix for counting all unique data points

for k=1:1:m; % using loop cycle for unique data points all combinations; for this k is used
    for l=1:1:m;
        for i=1:1:n-1;
          j=i+1;
            if x(i)== z(k) & x(j)==z(l);
              N(k,l) = N(k,l)+1;
            end
         i=i+1;
        end
     l=l+1;
    end
  k=k+1;
end
N
s=sum(N,2);
Tpm= N./s %transition probability matrix

Upvotes: 1

xvan
xvan

Reputation: 4855

%%Sample matrix
p=magic(8)
%%Fill rows and cols 3,5 with 0's
p([3 5],:)=0
p(:,[3 5])=0

%%The code

lb=[]
for k = [length(p):-1:1]
  if any(p(k,:)) | any(p(:,k))
     lb=[ [k],lb ]
  else
     p(k,:)=[]
     p(:,k)=[]
  end
end

lb keeps your original index

Upvotes: -1

Yvon
Yvon

Reputation: 2983

Step 1 - organize data and generate empty transition table

x= [12,14,6,15,15,15,15,6,8,8,18,18,14,14] 
xind = zeros(1,length(x));
u = unique(x) % find unique elements and sort
for ii = 1:length(u)
    xmask = x==u(ii); % locate all elements of a single value
    xind = xind+ii*xmask; % number them in the order listed in u
end

Output is labeled Markov chain (elements are labels instead of meaningful values)

>> u

u =

     6     8    12    14    15    18

>> xind

xind =

     3     4     1     5     5     5     5     1     2     2     6     6     4     4

Step 2 - build "from-to" table for each hop

>> T = [xind(1:end-1);xind(2:end)]

T =

     3     4     1     5     5     5     5     1     2     2     6     6     4
     4     1     5     5     5     5     1     2     2     6     6     4     4

Each column is a transition. First row is "from" label, second row is "to" label.

Step 3 - count frequencies and create transition table

p = zeros(length(u));
for ii = 1:size(T,2)
    px = T(1,ii); % from label
    py = T(2,ii); % to label
    p(px,py) = p(px,py)+1;
end

Output is aggregated frequency table. Each element is counts of a hop. Row number is "from" and column number is "to". 

>> p

p =

     0     1     0     0     1     0
     0     1     0     0     0     1
     0     0     0     2     0     0
     2     0     0     1     0     0
     1     0     0     0     3     0
     0     0     0     1     0     1

For example the 3 means 3 transitions from 5th label to 5th label (actual value is 15 to 15)

Step 4 - normalize row vectors to get probability table

>> p./repmat(sum(p,2),1,length(u))

ans =

         0    0.5000         0         0    0.5000         0
         0    0.5000         0         0         0    0.5000
         0         0         0    1.0000         0         0
    0.5000         0         0    0.5000         0         0
    0.2500         0         0         0    0.7500         0
         0         0         0    0.5000         0    0.5000

alternative loop version

for ii = 1:size(p,1)
    count = sum(p(ii,:));
    p(ii,:) = p(ii,:)/count;
end

Upvotes: 1

Related Questions