function [if_is if_mc po_is po_mc pbs] = bs(data,y0,covs,chain,L,sk)

par=get_pars(chain,sk);
is=impsamp(par,L,sk);

if_is=eval_imp(is,par,sk);
if_mc=eval_imp(chain,par,sk);
po_is=eval_pos(data,y0,covs,is,sk);
po_mc=eval_pos(data,y0,covs,chain,sk);


% M=size(chain.al,1);

ouli=find(po_mc.lpri<-1e20);
po_mc.lik(ouli)=[];
po_mc.llik(ouli)=[];
po_mc.pri(ouli)=[];
po_mc.lpri(ouli)=[];
if_mc.act(ouli)=[];
if_mc.log(ouli)=[];
M=length(po_mc.lik);

p0=mean((po_is.lik.*po_is.pri)./if_is.act);
disp(['Starting point ',num2str(p0,'%10.5e\n')])
pbs(1)=p0;

para=1e-3; dis=1; ite=1; tol=1e3;

while dis>=para && ite<=tol
    ite=ite+1;
    nume = mean(exp(po_is.llik + po_is.lpri)./(L*if_is.act + M*exp(po_is.llik + po_is.lpri)/pbs(ite-1)));
    deno = mean(if_mc.act ./(L*if_mc.act + M*exp(po_mc.llik + po_mc.lpri)/pbs(ite-1)));
    pbs(ite)=nume/deno;
    dis=abs((pbs(ite)-pbs(ite-1))/pbs(ite-1));
end
ite
disp(pbs(length(pbs)))