function [draws] = f_sample_sel(X1,X2,res,R,atmean,sampler)
    
    % infer dimensions
    H = size(res.beta_kernel,1);
    n = size(X1,1);
    
    [beta_mixing_sim,beta_kernel_sim,tau_sim] = f_simul_VB(res,R);
     beta_kernel_sim = permute(beta_kernel_sim,[2,3,1]);%repmat(res.beta_kernel,1,1,R);%permute(beta_kernel_sim,[2,3,1]);
     beta_mixing_sim = permute(beta_mixing_sim,[2,3,1]); %repmat(res.beta_mixing,1,1,R);%permute(beta_mixing_sim,[2,3,1]);
     tau_sim = tau_sim; %repmat(res.tau,R,1);%tau_sim;

     % repmat(res.beta_kernel,1,1,R);
     % repmat(res.beta_mixing,1,1,R);
     % repmat(res.tau,R,1);
    
    % initialize variables
    pred   = zeros(n,1);
    pii    = zeros(H,R);
    nu     = zeros(H,1);
    cum_nu = zeros(H,1);
    draws  = zeros(R,1);
    
    % cycle through the observations
    for r = 1:R
        if atmean == 0
            beta_mixing = (beta_mixing_sim(:,:,r));
            beta_kernel = (beta_kernel_sim(:,:,r));
            tau         = tau_sim(r,:);
        elseif atmean == 1
            beta_mixing = mean(beta_mixing_sim,3);
            beta_kernel = mean(beta_kernel_sim,3);
            tau         = mean(tau_sim);    
        elseif atmean == 2
            beta_mixing = res.beta_mixing;
            beta_kernel = res.beta_kernel;
            tau         = res.tau;  
        end
        % first mixture component
        nu(1)     = 1/(1+exp(-(X2*beta_mixing(1,:)')));
        pii(1,r)    = nu(1);
        cum_nu(1) = 1-nu(1);
       
       
        % update other mixture components until H-1
        for h = 2:H-1
            nu(h)     = 1/(1+exp(-(X2*beta_mixing(h,:)')));
            pii(h,r)    = nu(h) * cum_nu(h-1);
            cum_nu(h) = (1-nu(h))*cum_nu(h-1);
        end
    
         % update the last mixture component
        nu(H)  = 1;

        if sampler == 1
            pii(H,r) = nu(H)*cum_nu(H-1);
            pii(:,r) = pii(:,r).*(round(tau,3)'~=0);
            pii(:,r) = pii(:,r)./sum(pii(:,r));
        elseif sampler == 3
            pii(H,r) = nu(H)*cum_nu(H-1);
            pii(:,r) = pii(:,r).*(round(tau,2)'~=0);
            pii(:,r) = pii(:,r)./sum(pii(:,r));
        elseif sampler == 2
            pii(H,r) = nu(H)*cum_nu(H-1);
            if any(pii(:,r)<.05)
                pii((pii(:,r))<.05,r) = 0;
                pii(:,r) = pii(:,r)./sum(pii(:,r));
            end
            pii(:,r) = pii(:,r)./sum(pii(:,r));
        end


        G(r) = find(mnrnd(1,pii(:,r)));

        draws(r) = X1*beta_kernel(G(r),:)' + sqrt(1/(tau(G(r))))*randn;


    end
    
   
end




