function [v_mri,v_si,v_a_mri,v_a_si] = bayes_nonpar_mri_si(theta,alpha,ta);

% function [v_mri,v_a_mri,v_si,v_a_si] = bayes_nonpar_mri_si(theta,alpha,ta);

sz = size(theta);
nsub = 8;%sz(end);
nsamp = sz(1);%length(theta);
ntest = 4;
baserate = 1;%[1 2 3 4 5];

% theta
if ta(1)

    for br = baserate
        for subj = 1:nsub
            for mri = 1:ntest
                disp([br subj mri])
                for samp = 1:nsamp
                    if mri==1
                        pa_aa = theta(samp,1,1,br,subj)+theta(samp,2,1,br,subj);
                        pa_ab = theta(samp,1,2,br,subj)+theta(samp,2,2,br,subj);
                        dif = pa_aa - pa_ab;
                        v_mri(samp,mri,3,br,subj) = pa_aa;
                        v_mri(samp,mri,2,br,subj) = pa_ab;
                        v_mri(samp,mri,1,br,subj) = dif;
                    elseif mri==2
                        pb_ba = theta(samp,3,3,br,subj)+theta(samp,4,3,br,subj);
                        pb_bb = theta(samp,3,4,br,subj)+theta(samp,4,4,br,subj);
                        dif = pb_ba - pb_bb;
                        v_mri(samp,mri,3,br,subj) = pb_ba;
                        v_mri(samp,mri,2,br,subj) = pb_bb;
                        v_mri(samp,mri,1,br,subj) = dif;
                    elseif mri==3
                        p_aaa = theta(samp,1,1,br,subj)+theta(samp,3,1,br,subj);
                        p_aba = theta(samp,1,3,br,subj)+theta(samp,3,3,br,subj);
                        dif = p_aaa - p_aba;
                        v_mri(samp,mri,3,br,subj) = p_aaa;
                        v_mri(samp,mri,2,br,subj) = p_aba;
                        v_mri(samp,mri,1,br,subj) = dif;
                    elseif mri==4
                        p_bab = theta(samp,2,2,br,subj)+theta(samp,4,2,br,subj);
                        p_bbb = theta(samp,2,4,br,subj)+theta(samp,4,4,br,subj);
                        dif = p_bab - p_bbb;
                        v_mri(samp,mri,3,br,subj) = p_bab;
                        v_mri(samp,mri,2,br,subj) = p_bbb;
                        v_mri(samp,mri,1,br,subj) = dif;
                    end
                end
            end
        end
    end


    for br = baserate
        for subj = 1:nsub
            for si = 1:ntest
                disp([br subj si])
                for samp = 1:nsamp
                    if si==1
                        paaaa = theta(samp,1,1,br,subj);
                        pa_aa = theta(samp,1,1,br,subj)+theta(samp,2,1,br,subj);
                        p_aaa = theta(samp,1,1,br,subj)+theta(samp,3,1,br,subj);
                        v_si(samp,si,3,br,subj) = paaaa;
                        v_si(samp,si,2,br,subj) = pa_aa*p_aaa;
                        v_si(samp,si,1,br,subj) = paaaa-(pa_aa*p_aaa);
                    elseif si==2
                        pabab = theta(samp,2,2,br,subj);
                        pa_ab = theta(samp,1,2,br,subj)+theta(samp,2,2,br,subj);
                        p_bab = theta(samp,2,2,br,subj)+theta(samp,4,2,br,subj);
                        v_si(samp,si,3,br,subj) = pabab;
                        v_si(samp,si,2,br,subj) = pa_ab*p_bab;
                        v_si(samp,si,1,br,subj) = pabab-(pa_ab*p_bab);
                    elseif si==3
                        pbaba = theta(samp,3,3,br,subj);
                        pb_ba = theta(samp,3,3,br,subj)+theta(samp,4,3,br,subj);
                        p_aba = theta(samp,1,3,br,subj)+theta(samp,3,3,br,subj);
                        v_si(samp,si,3,br,subj) = pbaba;
                        v_si(samp,si,2,br,subj) = pb_ba*p_aba;
                        v_si(samp,si,1,br,subj) = pbaba-(pb_ba*p_aba);
                    elseif si==4
                        pbbbb = theta(samp,4,4,br,subj);
                        pb_bb = theta(samp,4,4,br,subj)+theta(samp,3,4,br,subj);
                        p_bbb = theta(samp,2,4,br,subj)+theta(samp,4,4,br,subj);
                        v_si(samp,si,3,br,subj) = pbbbb;
                        v_si(samp,si,2,br,subj) = pb_bb*p_bbb;
                        v_si(samp,si,1,br,subj) = pbbbb-(pb_bb*p_bbb);
                    end
                end
            end
        end
    end
else
   v_mri = 0;
   v_si = 0;
end

if ~ta(2)
    v_a_mri = 0;
    v_a_si = 0;
else
    % alpha
    for samp = 1:nsamp
        for br = baserate
            for i = 1:4
                alpha_p(samp,:,i,br) = alpha(samp,:,i,br)/sum(alpha(samp,:,i,br));
            end
            for mri = 1:ntest
                if mri==1
                    pa_aa = alpha_p(samp,1,1,br)+alpha_p(samp,2,1,br);
                    pa_ab = alpha_p(samp,1,2,br)+alpha_p(samp,2,2,br);
                    dif = pa_aa - pa_ab;
                    v_a_mri(samp,mri,3,br) = pa_aa;
                    v_a_mri(samp,mri,2,br) = pa_ab;
                    v_a_mri(samp,mri,1,br) = dif;
                elseif mri==2
                    pb_ba = alpha_p(samp,3,3,br)+alpha_p(samp,4,3,br);
                    pb_bb = alpha_p(samp,3,4,br)+alpha_p(samp,4,4,br);
                    dif = pb_ba - pb_bb;
                    v_a_mri(samp,mri,3,br) = pb_ba;
                    v_a_mri(samp,mri,2,br) = pb_bb;
                    v_a_mri(samp,mri,1,br) = dif;
                elseif mri==3
                    p_aaa = alpha_p(samp,1,1,br)+alpha_p(samp,3,1,br);
                    p_aba = alpha_p(samp,1,3,br)+alpha_p(samp,3,3,br);
                    dif = p_aaa - p_aba;
                    v_a_mri(samp,mri,3,br) = p_aaa;
                    v_a_mri(samp,mri,2,br) = p_aba;
                    v_a_mri(samp,mri,1,br) = dif;
                elseif mri==4
                    p_bab = alpha_p(samp,2,2,br)+alpha_p(samp,4,2,br);
                    p_bbb = alpha_p(samp,2,4,br)+alpha_p(samp,4,4,br);
                    dif = p_bab - p_bbb;
                    v_a_mri(samp,mri,3,br) = p_bab;
                    v_a_mri(samp,mri,2,br) = p_bbb;
                    v_a_mri(samp,mri,1,br) = dif;
                end
            end
        end
    end

    for samp = 1:nsamp
        for br = baserate
            for si = 1:ntest
                if si==1
                    paaaa = alpha_p(samp,1,1,br);
                    pa_aa = alpha_p(samp,1,1,br)+alpha_p(samp,2,1,br);
                    p_aaa = alpha_p(samp,1,1,br)+alpha_p(samp,3,1,br);
                    v_a_si(samp,si,3,br) = paaaa;
                    v_a_si(samp,si,2,br) = pa_aa*p_aaa;
                    v_a_si(samp,si,1,br) = paaaa-(pa_aa*p_aaa);
                elseif si==2
                    pabab = alpha_p(samp,2,2,br);
                    pa_ab = alpha_p(samp,1,2,br)+alpha_p(samp,2,2,br);
                    p_bab = alpha_p(samp,2,2,br)+alpha_p(samp,4,2,br);
                    v_a_si(samp,si,3,br) = pabab;
                    v_a_si(samp,si,2,br) = pa_ab*p_bab;
                    v_a_si(samp,si,1,br) = pabab-(pa_ab*p_bab);
                elseif si==3
                    pbaba = alpha_p(samp,3,3,br);
                    pb_ba = alpha_p(samp,3,3,br)+alpha_p(samp,4,3,br);
                    p_aba = alpha_p(samp,1,3,br)+alpha_p(samp,3,3,br);
                    v_a_si(samp,si,3,br) = pbaba;
                    v_a_si(samp,si,2,br) = pb_ba*p_aba;
                    v_a_si(samp,si,1,br) = pbaba-(pb_ba*p_aba);
                elseif si==4
                    pbbbb = alpha_p(samp,4,4,br);
                    pb_bb = alpha_p(samp,4,4,br)+alpha_p(samp,3,4,br);
                    p_bbb = alpha_p(samp,2,4,br)+alpha_p(samp,4,4,br);
                    v_a_si(samp,si,3,br) = pbbbb;
                    v_a_si(samp,si,2,br) = pb_bb*p_bbb;
                    v_a_si(samp,si,1,br) = pbbbb-(pb_bb*p_bbb);
                end
            end
        end
    end
end