function [state] = sample_chain(phi, Phi, T)

state(1) = sample_multi(phi);

for t=2:T
    
    state(t) = sample_multi(Phi(:,state(t-1)));
    
end

