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;