n = 1000000; observe = [1 1 1 0 0 1 1 1 0 1]; % Naive way % 1 state = floor(rand(n,1) * 3); while size(state,2) < length(observe) state = [state, mod(state(:,end) - (rand(n,1) < .3), 3)]; end obs = rand(size(state)) < state * .4 + .1; % 2 antecedent = ones(n,1); for i = 1 : length(observe) antecedent = antecedent & obs(:,i) == observe(i); end % Less naive way state = floor(rand(n,1) * 3); while size(state,2) < length(observe) state = [state, mod(state(:,end) - (rand(n,1) < .3), 3)]; end weight = ones(n,1); for i = 1 : length(observe) weight = weight .* ((1 - observe(i)) + ... (state(:,i) * .4 + .1) * (observe(i) * 2 - 1)); end % Even less naive way, forgetting all but the last state so we % actually have room to enumerate all states exhaustively state = (0:2)'; weight = repmat(1/3, 3, 1) .* ((1 - observe(1)) + ... (state * .4 + .1) * (observe(1) * 2 - 1)); for i = 2 : length(observe) state = [state; mod(state - 1, 3)]; weight = [weight * .7; weight * .3]; weight = full(sparse(state + 1, ones(size(state)), weight)); state = (0 : length(weight)-1)'; weight = weight .* ((1 - observe(i)) + ... (state * .4 + .1) * (observe(i) * 2 - 1)); end