function [plower, pupper] = findbounds(NR, n, jconsec)
%Matlab code for "A Set Probability Technique for Detecting
%Relative Time Order Across Multiple Neurons", Neural Computation,
%In Press.
%Takes 3 arguments (NR, n, jconsec)
% NR is the reference sequence length
% n is the word length
% j is the number of consecutive increasing elements
%Returns plower and pupper = upper and lower bounds on
%Pr(j or more strictly increasing)
%If an exact solution is possible then pupper=plower
%Anne C Smith April, 2005
res = [];
cnt = 0;
%compute pj = Pr(j in a row or more are increasing)
j = jconsec;
%Eq. 1
pj = factorial(NR)/factorial(j)/factorial(NR-j)/(NR^j);
%3 cases depending on relative sizes of n and j
if(n==j)
pfinal = pj; %subsequence length same as sequence
gfinal = pj;
elseif( 2*j>n ) %case of a very short sequence-exact
pjp1 = factorial(NR)/factorial(j+1)/factorial(NR-j-1)/(NR^(j+1));
pfinal = pj + (n-j)*(pj-pjp1); % Eq. 8
gfinal = pfinal; % upper and lower bound the same
elseif( 2*j<=n ) %deals with longer sequences
pjp1 = factorial(NR)/factorial(j+1)/factorial(NR-j-1)/(NR^(j+1));
f = []; g = []; % f is upper bound, g is lower bound
f(1) = pj;
g(1) = pj;
for r = 2:j
f(r) = pj-pjp1;
g(r) = f(r);
end
for r = j+1:n-j+1;
xx = 1:r-j-1;
s = sum(f(xx));
f(r) = (pj-pjp1) * (1-s); %upper bound, Eq. 15
g(r) = (pj-pjp1) - (r-j)*pj*pj; %lower bound, Eq. 16
if(g(r) < 0)
g(r) = 0;
end
end
pfinal = sum(f);
gfinal = sum(g);
end
if(pfinal > 1) % check pfinal not slightly over 1
pupper = 1;
else
pupper = pfinal;
end
plower = gfinal;