INTELLIGENT WORK FORUMS FOR ENGINEERING PROFESSIONALS
Come Join Us!
Are you an Engineering professional? Join EngTips now!
 Talk With Other Members
 Be Notified Of Responses
To Your Posts
 Keyword Search
 OneClick Access To Your
Favorite Forums
 Automated Signatures
On Your Posts
 Best Of All, It's Free!
*EngTips's functionality depends on members receiving email. By joining you are opting in to receive email.
Posting Guidelines
Promoting, selling, recruiting, coursework and thesis posting is forbidden.

Statetransition matrix(3)

BrCo (Marine/Ocean) (OP) 
2 Feb 09 18:06 
Dear All, I would appreciate your help. I have a vector (e.g., x=[4 6 2 3 3 1 5 3 4] ) I need to calculate the statetransition matrix, that is, how many times does each number follow each number. In other words, what are the chances of 3 to come after 1, what are the chances of 2 to come after 5, etc. This may seem easy in a small vector, but my vectors each have about 20,000 numbers, so an efficient system will help a lot. Thanks! 

A brute force program in Octave for 20000 numbers 16 took 2.2 seconds to run and 2 minutes to write. The only line of interest in the entire exercise is ++x(z(i),z(i+1)); and that is about as interesting as lettuce. Cheers
Greg Locock
SIG:Please see FAQ731376: EngTips.com Forum Policies for tips on how to make the best use of EngTips. 

The same prog in matlab runs in 0.016 s, although I had to rewrite the interesting line because matlab doesn't like ++. So, how fast, exactly, do you need this algorithm to be? Cheers
Greg Locock
SIG:Please see FAQ731376: EngTips.com Forum Policies for tips on how to make the best use of EngTips. 

FeX32 (Mechanical) 
2 Feb 09 23:37 
Histogram it is. My first thought was this crude code: clear clc x=1:6; for i=1:2000 y(i)=randsample(x,1); end su=2000/6; z1=zeros(1,6); z2=z1;z3=z2;z4=z3;z5=z4;z6=z5; for i=1:1999 if y(i)==1 for k=1:6 if y(i+1)==k z1(k)=z1(k)+1; end end end if y(i)==2 for k=1:6 if y(i+1)==k z2(k)=z2(k)+1; end end end if y(i)==3 for k=1:6 if y(i+1)==k z3(k)=z3(k)+1; end end end if y(i)==4 for k=1:6 if y(i+1)==k z4(k)=z4(k)+1; end end end if y(i)==5 for k=1:6 if y(i+1)==k z5(k)=z5(k)+1; end end end if y(i)==6 for k=1:6 if y(i+1)==k z6(k)=z6(k)+1; end end end end figure(1) bar(x,z1), title('followers of 1') figure(2) bar(x,z2), title('followers of 2') figure(3) bar(x,z3), title('followers of 3') figure(4) bar(x,z4), title('followers of 4') figure(5) bar(x,z5), title('followers of 5') figure(6) bar(x,z6), title('followers of 6') It's quick and was simple to write. From each z (ie. z1 to z6) you have the number of times each number follows it. Knowing the total population size and the variation in numbers you can easily get a percentage of what number follows what. Greg I'm sure your code is better. But I'm not sure how you utilized that interesting line. Fe 

FeX32 (Mechanical) 
2 Feb 09 23:41 
Also should be able to put those 6 'if' loops in a single 'for' loop to compress the size. just a thought, Fe 

num_spots=input('How many spots '); num_trials=input('How many trials '); for i=1:num_spots for j=1:num_spots x(i,j)=0; end end for i=1:num_trials z(i)=1+floor(num_spots*rand(1)); end tic for i=1:(num_trials1) x(z(i),z(i+1))=x(z(i),z(i+1))+1;%was the interesting line end toc x Cheers
Greg Locock
SIG:Please see FAQ731376: EngTips.com Forum Policies for tips on how to make the best use of EngTips. 

BrCo (Marine/Ocean) (OP) 
3 Feb 09 1:20 

FeX32 (Mechanical) 
3 Feb 09 17:55 
Glad to help. By the way, I always thought the term 'statetransition matrix' was used in control engineering. Fe 

That's why I didn't reply. Had I known it was a regular ML golf problem, I may have played.  Steve 

BrCo: This exact question came up on CSSM last year. There are many ways to solve these sorts of problems. My favourite is to use the sparse() function: % Some data with repeated sequences x=[1 6 1 6 4 4 4 3 1 2 2 3 4 5 4 5 2 6 2 6 2 6]; sparse(x(1:end1),x(2:end),1) ans = (3,1) 1 (6,1) 1 (1,2) 1 (2,2) 1 (5,2) 1 (6,2) 2 (2,3) 1 (4,3) 1 (3,4) 1 (4,4) 2 (5,4) 1 (6,4) 1 (4,5) 2 (1,6) 2 (2,6) 3 Or if you want it in array form, where A(i,j) holds the number of changes from i to j: A=full(sparse(x(1:end1),x(2:end),1)) A = 0 1 0 0 0 2 0 1 1 0 0 3 1 0 0 1 0 0 0 0 1 2 2 0 0 1 0 1 0 0 1 2 0 1 0 0  Steve 

FeX32 (Mechanical) 
4 Feb 09 11:22 
Nice and simple. I never used sparse before. What is CSSM? Fe 

It was my answer in CSSM, so I'm not plagiarising, just copying.  Steve 

Sorry, forgot a star for you. Slow isn't it? where is that construction of x(1:end1),x(2:end) documented? What I mean is how does it 'know' that it can and should increment the pointer in each of the two x() statements together? Cheers
Greg Locock
SIG:Please see FAQ731376: EngTips.com Forum Policies for tips on how to make the best use of EngTips. 

Greg, Just look through the documentation for sparse(). In the form I've used, x(1:end1) and x(2:end) are just two vector inputs representing the row and column for each of the elements in the third vector. For example, this produces a sparse version of eye(3): sparse([1 2 3],[1 2 3],[1 1 1]) If the values of the elements are all the same, a single value can be (re)used: sparse([1 2 3],[1 2 3],1) This sort of assignment can't be done with normal arrays, there is no equivalent syntax (you run into all that sub2ind, ind2sub confusion). Now the trick: The neat (documented, but not obvious) thing about sparse()is that if a coordinate pair is repeated, the values assigned to it are summed, so this: sparse([3 1 2 3 3],[3 1 2 3 3],1) gives (1,1) = 1 (2,2) = 1 (3,3) = 3 So x(1:end1) are the "from" coordinates and x(2:end) are the "to" coordinates. Each (from,to) pair is assigned a value of 1. Each repeat of a (from,to) pair increments the value at that location. This is what I love about Matlab, sometimes it all just falls together.  Steve 

FeX32 (Mechanical) 
5 Feb 09 9:03 
Thanks for the info on CSSM. I wish I knew sparse existed before TTFN Fe 

Ah, ok it is a 'sparse' thing, not an indexing mode that I have yet to get to grips with (ie almost all of them). Cheers
Greg Locock
SIG:Please see FAQ731376: EngTips.com Forum Policies for tips on how to make the best use of EngTips. 



