Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
calculating markov transition probability

Subject: calculating markov transition probability

From: Yidian

Date: 24 Apr, 2013 10:33:07

Message: 1 of 6

I am calculating markov transition probability using matlab.
basically, i am solving A*x=B, where A is [1 x n], x is [n x n], B is [1 x n].
the constraints are
1. x is an upper triangle matrix, i.e. [x11, x12, x13; 0 x22, x23; 0 0 x33]
2. each row in x sums to 1
3. x33 ==1

does anyone know a way to solve it? any input would be helpful!

Subject: calculating markov transition probability

From: Bruno Luong

Date: 24 Apr, 2013 13:45:09

Message: 2 of 6

% Fake data
n = 3;
X = triu(rand(n));
X = bsxfun(@rdivide,X,sum(X,2))
A = rand(1,n);
B = A*X;
clear X

% Engine
M = kron(eye(n),A);
% Matrix to constraint lower part is zero, C1*X(:) = b1
iz = find(tril(ones(n),-1));
nz = length(iz);
n2 = size(M,2);
C1 = accumarray([(1:nz)' iz], 1, [nz n2]);
b1 = zeros(nz,1);
% Matrix to constraint column-sums X are 1
C2 = kron(ones(1,n),eye(n));
b2 = ones(n,1);
% Merge all the constrains
C = [C1; C2];
b = [b1; b2];
% Resolution A*X = B, with constraint C*X(:) = b
K = [M'*M C'; C zeros(size(C,1))];
RHS = [M'*B(:); b];
Y = pinv(K) * RHS;
Y = Y(1:n^2);
Y = reshape(Y,[n n]);
Y % should be X

% Bruno

Subject: calculating markov transition probability

From: Bruno Luong

Date: 24 Apr, 2013 16:03:09

Message: 3 of 6

Another way of resolution is (the rest is the same):

% Resolution
Q = null(C);
Y0 = C\b;
Y = Y0 + Q*(pinv(M*Q) * (B(:)-M*Y0));
Y = reshape(Y,[n n])

% Bruno

Subject: calculating markov transition probability

From: Matt J

Date: 24 Apr, 2013 16:18:07

Message: 4 of 6

"Yidian" wrote in message <kl8cd3$q3h$1@newscl01ah.mathworks.com>...
> I am calculating markov transition probability using matlab.
> basically, i am solving A*x=B, where A is [1 x n], x is [n x n], B is [1 x n].
> the constraints are
> 1. x is an upper triangle matrix, i.e. [x11, x12, x13; 0 x22, x23; 0 0 x33]
> 2. each row in x sums to 1
> 3. x33 ==1


What about xij>=0? No constraint like that?

Subject: calculating markov transition probability

From: Bruno Luong

Date: 24 Apr, 2013 18:23:09

Message: 5 of 6

"Yidian" wrote in message <kl8cd3$q3h$1@newscl01ah.mathworks.com>...
> I am calculating markov transition probability using matlab.
> basically, i am solving A*x=B, where A is [1 x n], x is [n x n], B is [1 x n].
> the constraints are
> 1. x is an upper triangle matrix, i.e. [x11, x12, x13; 0 x22, x23; 0 0 x33]
> 2. each row in x sums to 1
> 3. x33 ==1
>
> does anyone know a way to solve it? any input would be helpful!

Note that this problem is undetermined.

If X one solution, { X+s*P } are set of solution, where P is (3x3) defined from A as

P := [0 A(2) -A(2);
        0 -A(1) A(1);
        0 0 0].

In fact it is easy to show that:
1. (A*P) = 0,
2. sum(P,2) = 0, and
3. P is upper triangular.

% Bruno

Subject: calculating markov transition probability

From: Bruno Luong

Date: 25 Apr, 2013 07:01:10

Message: 6 of 6

With the kernel characterized, we can use "\" instead of pinv.

% Engine
M = kron(eye(n),A);
% Matrix to constraint lower part is zero, C1*X(:) = b1
iz = find(tril(ones(n),-1));
nz = length(iz);
n2 = size(M,2);
C1 = accumarray([(1:nz)' iz], 1, [nz n2]);
b1 = zeros(nz,1);
% Matrix to constraint column-sums X are 1
C2 = kron(ones(1,n),eye(n));
b2 = ones(n,1);
% Merge all the constrains
C = [C1; C2];
b = [b1; b2];
% kernel of X -> A*X
P = [0 A(2) -A(2);
     0 -A(1) A(1);
     0 0 0];
Q = null([C; P(:)']);
X0 = C\b;
X = X0 + Q*((M*Q) \ (B(:)-M*X0));
X = reshape(X,[n n]);
% Generate a random solution
s = rand();
X = X + s*P;
norm(A*X-B)

% Bruno

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us