get_slope_ests.m

function ests = get_slope_ests(N, filtlen)
%Missing values are dealt with by sampling from the Nratios
%N is a row vector of raw data (not log transformed)
%Following Holmes 2001 but mu estimated by means  of runsums not the
%regression

%get rid of missing values and zeros
if(ismember(0,N)) N=N+1; end
Nratios = N(2:end)./N(1:end-1);
Nratios = Nratios(Nratios > 0);
for(j = 2:length(N))
    if(N(j)==-99) N(j)= N(j-1)*sample(1,Nratios); end %interp by sampling from Nt+1/Nt ratios; sample is a subfunction below
end

Nratios = log(N(2:end)./N(1:(end-1)));

switch filtlen
    case 4,
        runsum = N(1:(end-3))+N(2:(end-2))+N(3:(end-1))+N(4:end);
    case 3,
        runsum = N(1:(end-2))+N(2:(end-1))+N(3:end);
    case 2,
        runsum = N(1:(end-1))+N(2:(end));
    case 1,
        runsum = N(1:end);
end
runlen = length(runsum);

Rratios = log(runsum(2:end)./runsum(1:(end-1)));
mu = mean(Rratios(1:end));
%Get sigma2p using slope of variance versus tau
%This uses var(tmp) which is sensitive to outliers, but the mad(x) is inefficient and data are
%often limited.  But if normality problems are evident and data are sufficient, try (1.3*mad(tmp))^2 as a variance estimator
taumax=3;  %4 works fine too
varrunsums=[];
for j=1:taumax,
    tmp = log(runsum((j+1):runlen))-log(runsum(1:(runlen-j)));
    varrunsums = [varrunsums var(tmp)];
end
%Get variance slope by linear regression
%alternatively you can use just the first and last points: sigma2slp=(varrunsums(4)-varrunsums(1))/3
[a b c d stats]=regress(varrunsums',[ones(taumax,1) (1:taumax)']);
s2=max(10^(-5),a(2));  %don't allow negative values
s2np = max(0,(var(Nratios)-s2)/2);

ests.s2p=s2;
ests.mu=mu;
ests.s2s=s2np;

function samp=sample(N,X)
%sample without replacement
samp=zeros(1,N);
I=zeros(1,N);
for(i=1:N)
    newI=unidrnd(length(X));
    while(ismember(newI,I)) newI=unidrnd(length(X)); end
    I(i)=newI;
    samp(i)=X(newI);
end