function [pL ,pU]=gen(a,b,t,int)
%function [pL ,pU]=gen(a,b,t,int)
%a: vector of lower bounds
%b: vector of upper bounds
%t: vector of information times usually 1,2,3...k
%int: vector of number of intervals to use for the integration
%OUTPUT
%pL=Cumulative probability of y(k) going below a(k)
%pU=Cumulative probability of y(k) going above b(k)
%y is standard normal under Ho
%y=z(t)/sqrt(t) where Cov(z(t_1),z(t_2))=min(t_1,t_2);
%example null hypothesis, Ware and Demets Biometrika 69(3) pg 663 table line 9.
%[pl,pu]=gen(-1.584./(1.651*sqrt(1:3))+.5*1.651*sqrt(1:3),1.988*sqrt(3./(1:3)),1:3,500*ones(1,3))
%example alternative: subtract mean from a and b, note that mean increases with t
%[pl,pu]=gen(-1.584./(1.651*sqrt(1:3))+.5*1.651*sqrt(1:3)-sqrt(1:3)*2,1.988*sqrt(3./(1:3))-sqrt(1:3)*2,1:3,500*ones(1,3))
if (length(a)~=length(b))|length(a)~=length(t)|length(a)~=length(int) 
'ERROR WRONG MATRIX LENGTH'
end
d=(b-a)./int;m=length(a);sq2pi=sqrt(2*pi);H=1:int(1);E=ones(1,int(1));
xo=a(1)+((1:int(1))-.5*E)*d(1);
pU(1)=normcdf(-(sqrt(t(1))*b(1))/sqrt(t(1)));
M=((d(1)/sq2pi)*exp(-(sqrt(t(1))*xo).^2/(2*t(1))))';
pL(1)=normcdf(sqrt(t(1))*a(1)/sqrt(t(1)));
for k=2:m
   VU=normcdf(-(sqrt(t(k))*b(k)*E-sqrt(t(k-1))*xo)/sqrt(t(k)-t(k-1)));
   VL=normcdf((sqrt(t(k))*a(k)*E-sqrt(t(k-1))*xo)/sqrt(t(k)-t(k-1)));
   pL(k)=pL(k-1)+VL*M;
   pU(k)=pU(k-1)+VU*M;
   x=a(k)+((1:int(k))-.5*ones(1,int(k)))*d(k);
   M=(d(k)*sqrt(t(k))/(sq2pi*sqrt(t(k)-t(k-1))))*...
       exp(-(sqrt(t(k))*(x'*ones(1,int(k-1)))-sqrt(t(k-1))*...
      (ones(int(k),1)*xo)).^2/(2*(t(k)-t(k-1))))*M;
       xo=x;
  end
%save c:\aa\matlab\gen.mat