% Basic Matlab code for forming the stopover capture-recapture likelihood

function L = stopover(theta,x)

% x: encounter histories
% D: number of observed individuals; T: number of encounter occasions

[D,T] = size(x);

for i = 1:D
    for t = 1:T
        if x(i,t)==1 && sum(x(i,1:t-1))==0
            f(i)=t;
        end
    end
    for t = T:-1:1
        if x(i,t)==1 && sum(x(i,t+1:T))==0
            l(i)=t;
        end
    end
end

% Parameters

% N: population size (super-population); n denotes the number of uncaught
% individuals

n = exp(theta(1));
N = n+D;

% beta(t): proportion of the N birds who are first available for capture at
% occasion t.  Note beta starts at beta(1) here rather than beta(0).

sumbeta = 0;
maxbeta = 4;
for t = 1:maxbeta-1
    sumbeta = sumbeta+exp(theta(t+1));
end

for t = 1:maxbeta-1
    beta(t) = exp(theta(t+1))/(1+sumbeta);
end

beta(maxbeta) = 1/(1+sumbeta);

for t = maxbeta+1:T
    beta(t) = 0;
end

% p(t,a): probability an individual is captured at occasion t given it
% arrived a occasions ago

for a = 1:T
    for t = 1:T
        p(t,a) = ilogit(theta(5));
    end
end

% phi(t,a): probability an individual in the study at occasion t remians in
% the study until occasion t+1, given that it arrived a occasions ago

for t = 1:T
    for a = 1:T
        phi(t,a) = ilogit(theta(t+5)+a*theta(15));
    end
end

% probabilities for observed individual i

for i = 1:D
    prob(i) = 0;
    for b = 1:f(i)
        for d = l(i):T
            prodphi1 = 1;
            for j = b:d-1
                prodphi1=prodphi1*phi(j,j-b+1);
            end
            prodp1 = 1;
            for j = b:d
                prodp1 = prodp1*p(j,j-b+1)^(x(i,j))*(1-p(j,j-b+1))^(1-x(i,j));
            end
            prob(i) = prob(i)+beta(b)*prodphi1*(1-phi(d,d-b+1))*prodp1;
        end
    end
end

% probability for unobserved individuals

prob0 = 0;
for b = 1:T
    for d = b:T
        prodphi2 = 1;
        for j = b:d-1
            prodphi2=prodphi2*phi(j,j-b+1);
        end
        prodp2 = 1;
        for j = b:d
            prodp2 = prodp2*(1-p(j,j-b+1));
        end
        prob0 = prob0+beta(b)*prodphi2*(1-phi(d,d-b+1))*prodp2;
    end
end

prodprob=0;
for i = 1:D
    prodprob = prodprob + log(prob(i));
end

L = gammaln(N+1) - gammaln(n+1) + prodprob + n*log(prob0);

L = -L;