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
