function site = BootRecycling(site,Psigb,M,mamflag)
%Calculates the pairs of [theta L(theta)]
%The output is stored in the same structured array sent as input
%from eq. 4.1 of Davison et al., 1995

N=3*M;
for jj=1:length(site)

start=site(jj).start;
de=site(jj).de;
n=length(de(:,1));

b2Pmatrix=zeros(N,n);
b2mamvec=zeros(N,1);
pairs=zeros(M,2);

%make the kernel density function for p*(theta)
switch mamflag
    case 1
        mam=getmam3(de,Psigb(1),start);
    case 2
        mam=getmam3ul(de,Psigb(1),start);
end

h=mam(2); %set bandwidth to mam se

%build the database of N second-level bootstraps
for i=1:N
    sigb=randn*Psigb(2)+Psigb(1);  % randomise sigma_b
    R = unidrnd(n,n,1);            % n random integers with replacement
    
    %record frequencies of each n, used in calc of product term
    for k=1:n
        b2Pmatrix(i,k)=length(find(R==k));
    end
    
    Rde=de(R,1:2);                 % select corresponding de values
    r = randn(n,1);              
    de2 = Rde(:,1);
    de2(:,2)=Rde(:,2); 

    switch mamflag
        case 1
            mam=getmam3(de2,sigb,start);
        case 2
            mam=getmam3ul(de2,sigb,start);
    end
    
    b2mamvec(i,1)=mam(1);
end

%Calculate the bootstrap likelihoods
%make a first-level bootstrap dataset
for b1=1:M
    b1Prowvec=[];
    sigb=randn*Psigb(2)+Psigb(1);  % randomise sigma_b
    R = unidrnd(n,n,1);            % n random integers with replacement
    
    Rde=de(R,1:2);                 % select corresponding de values
    r = randn(n,1);              
    de2 = Rde(:,1);
    de2(:,2)=Rde(:,2);
    


    switch mamflag
        case 1
            mam=getmam3(de2,sigb,start);
        case 2
            mam=getmam3ul(de2,sigb,start);
    end
    
    theta=mam(1);
    
    %make the kernel density function for p*(theta)
    w=ones(n,1);
    [f]=ksdensity(de2(:,1),de2(:,1),'width',h,'weights',w); %weighted kernel density
    pStarTheta=[f/sum(f)]';
    
    %calc likelihood------------------
    thetavec=repmat(theta,N,1);
    kdthis=(b2mamvec-thetavec)./h;
    kd1=normpdf(kdthis);
    
    %the product term:
    Pmat=pStarTheta/(1/n);
    Pmat=repmat(Pmat,N,1);
    prodterm=prod(Pmat.^b2Pmatrix,2);
    
    kd2=kd1.*prodterm;
    kd=sum(kd2);
    
    likelihood=(1/(N*h))*kd;
    pairs(b1,1:2)=[theta likelihood];
    site(jj).pairs=pairs;
end

end
