function [st_draw, smooth_t, filter_t]=state_draw_all_uni(y,x,beta1_draw,beta0_draw, sigma1_draw, sigma0_draw, tp11_draw, tp00_draw)

T=size(y,1);
st_draw=zeros(T,1);
U=unifrnd(0,1,T,1);

% first let's define some initial settings
P=[tp11_draw, 1-tp11_draw; 1-tp00_draw, tp00_draw]; % which is the transistion matrix
filter_t=zeros(T,2);       % which is the filter proability: filter_t(t,1)=Pr(St=1|Yt)
conditional_t=zeros(T,2);  % conditional_t=Pr(St=k|Yt-1)
smooth_t=zeros(T,2);       % smooth_t=Pr(St=k|Y_T)=Pr(St=k|Y_t)*Pr(St+1|St), which is the product of {filter probability}*{transtion probability}

filter_t(1,:)=[0.5, 0.5];  % assign the initial filter_t value
conditional_t(1,:)=[0.5, 0.5];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate the filter probability, then calculate  posterior conditional probability: P(st|Y_T)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:T
    likelihood(i,1)=1/sqrt(2*pi*sigma1_draw)*exp(-0.5*(y(i,1)-x(i,:)*beta1_draw).^2/(sigma1_draw));  % when st=1
    likelihood(i,2)=1/sqrt(2*pi*sigma0_draw)*exp(-0.5*(y(i,1)-x(i,:)*beta0_draw).^2/(sigma0_draw));  % when st=0
    
    if sum(likelihood(i,:)==0)>0        % i.e. two likelihood are both zero, avoid 0/0 problem
        likelihood(i,:)=1e-30.*ones(2,1);
    end
    conditional_t(i,:)=filter_t(i-1,:)*P;
    joint_likelihood(i,1)=sum(conditional_t(i,:).*likelihood(i,:));
    filter_t(i,:)=(conditional_t(i,:).*likelihood(i,:))./joint_likelihood(i,1);  % Normalized filter_t
end

st_draw(end)=(U(end)<=filter_t(end,1));  % if random uniform number is less than filter_t, we set it as 1, otherwise set it as 0
% caculate the posterior conditional probability: P(st|Y_T)
for i=T-1:-1:1
    if st_draw(i+1)==1
        smooth_t(i,:)=filter_t(i,:).*P(:,1)';                % -->[filter_t1*p11, filter_t0* p01]
    elseif st_draw(i+1)==0
        smooth_t(i,:)=filter_t(i,:).*P(:,2)';                % -->[filter_t1*p10, filter_t0* p00]
    end
    smooth_t(i,:)=smooth_t(i,:)./sum(smooth_t(i,:),2);   % Normalize smooth_t
    st_draw(i)=(U(i)<=smooth_t(i,1));
end

smooth_t(end,:)=filter_t(end,:);
end
