% post_sampling_f
% have sampled from the prior and have observed number of failures. Use the
% observations to give weight to each of the sampled values
for k=1:KK      %for each sampled value...
        
    l_weight(k)=0;
    s_r = samples(k,2) + samples(k,1)*exp(-samples(k,3).*(0:I));    %compute the failure rate for the sampled values 
    expected_fails = expected_function(samples(k,1),samples(k,2),samples(k,3));
    variance_fails = variances_func(samples(k,1),samples(k,2),samples(k,3));
    for i=2:I   %first obs is zero since components are working when they are produced.
        l_weight(k) = l_weight(k) + log_pnorm(obs(i),expected_fails(i),variance_fails(i));
    end
    
%    weight(k) =exp(l_weight);
end
clear s_r expected_fails variance_fails k
display('computing weights')

%maxt=max(exp(l_weight));
%temp1(1:KK) = l_weight(1:KK)-maxt;
%temp2 = exp(temp1);
%temp3= sum(temp2);
%temp4=log(temp3);
%l_s=maxt + temp4;

%weight(1:KK) = exp(l_weight(1:KK)-l_s)
%
maxt = max(l_weight);

for k=1:KK
    temp1(k) = l_weight(k) - maxt;
    temp2(k) = exp(temp1(k));
end
temp3 = sum(temp2);
temp4 = maxt+log(temp3);
s=exp(temp4);

for k=1:KK
    weight(k) = exp(l_weight(k))/s;
end
clear temp1 temp2 temp3 temp4 l_weight maxt
s_cdf=cumsum(weight);

for i=1:MM
    u=rand(1);  
    sample_x = 1;
        while (u > s_cdf(sample_x))
            sample_x=sample_x+1;
        end
    posterior_sim(i,:)=samples(sample_x,:);
end

posterior_sim_stored(count,:,:) = posterior_sim(:,:);
clear sample_x u i l_weight