c MAIN PROGRAM c This program reads in an integer-valued dataset, computes maximum likelihood c estimators for tau, theta0, and theta1, along with a 95% confidence interval c for tau. c double precision tau(200), theta1(200), theta0(200), lik(200), & taubest, t1best, t0best, alpha, llim, rlim, dT integer x(1000), T alpha=0.05d0 i=1 20 read(*,*,end=100) x(i) i=i+1 goto 20 100 T = i-1 call est1(T,x,theta0,theta1,tau,lik,taubest,t0best,t1best) call confint(100,x,T,alpha,llim,rlim) print*,'The maximum likelihood estimates of tau is',tau ptint*,'and a 95% confidence interval for tau is (',llim,rlim,')' print*,'The estimates of the Poisson rates are',t0best,'and',t1best, & 'before and after the change-point, respectively.' end c c A c Function to compute the first function of tau used for applying the quadratic c equation. c double precision function a(T,tau) double precision tau integer T double precision p a = (T-tau)*(T*p(tau)-tau) return end c c B c Function to compute the second function of tau used for applying the c quadratic equation. c double precision function b(T,tau,x) double precision tau integer T, x(T) double precision p integer brack double precision sum1, sum2 integer i c Compute sum of all X's sum1 = 0.d0 do 100 i = 1, T sum1 = sum1 + x(i) 100 continue c Compute sum of all X's _after_ the change-point sum2 = 0.d0 do 200 i = brack(tau)+2,T sum2 = sum2 + x(i) 200 continue b = -1*(T-tau)*p(tau)*sum1 - T*p(tau)*sum2 + tau*sum2 & + tau*(1.d0-p(tau))*x(brack(tau)+1) return end c c BRACK c Function to return the greatest integer function of tau. c integer function brack(tau) double precision tau brack = idint(tau) return end c c C c Function to compute the third function of tau used for applying the quadratic c equation. c double precision function c(T,tau,x) double precision tau integer T, x(T) double precision p integer brack double precision sum1, sum2 integer i c Compute sum of all X's sum1 = 0.d0 do 100 i = 1, T sum1 = sum1 + x(i) 100 continue c Compute sum of all X's _before_ the change-point sum2 = 0.d0 do 200 i = brack(tau)+2, T sum2 = sum2 + x(i) 200 continue c = p(tau)*sum1*sum2 return end c c CONFINT c Routine to compute a 100(1-alpha)% confidence interval for the change-point c tau. This works by computing all possible 100(1-alpha)% confidence intervals c (depending on the unit division variable unitdiv) and returns the one with c the shortest length. The final interval is (llimit,rlimit). c subroutine confint(unitdiv,x,T,alpha,llimit,rlimit) double precision alpha, llimit, rlimit integer unitdiv, x(T), T double precision xx, th0, th1, intgrl, baseline, & likfunc(100000), sofar integer ilikfunc, nlikfunc, ifirst, ilast, ifbest, ilbest double precision lnlik nlikfunc=(T-2)*unitdiv-1 c Compute the likelihood on an evenly-spaced grid intgrl = 0.d0 do 400 ilikfunc = 1,nlikfunc xx = 1.d0+(ilikfunc+0.d0)/(unitdiv+0.d0) call esttheta(T,xx,x,th0,th1) if(ilikfunc.eq.1) baseline = lnlik(T,xx,x,th0,th1) c Quick and dirty fix: if either estimate is zero, set the c likihood to zero. if(th0.eq.0.or.th1.eq.0) then likfunc(ilikfunc) = 0.d0 else likfunc(ilikfunc) = dexp(lnlik(T,xx,x,th0,th1)-baseline) endif intgrl = intgrl + likfunc(ilikfunc) 400 continue c Compute the first confidence interval ifirst = 1 sofar = likfunc(ifirst) ilast = 1 450 ilast = ilast + 1 sofar = sofar + likfunc(ilast) if(sofar/intgrl.lt.1.d0-alpha) goto 450 ifbest = ifirst ilbest = ilast c Compute the rest of the intervals, comparing lengths as it goes. 550 ifirst = ifirst + 1 sofar = sofar - likfunc(ifirst) if(sofar/intgrl.ge.1.d0-alpha) then if(ilast-ifirst.lt.ilbest-ifbest) then ifbest = ifirst ilbest = ilast endif goto 550 endif 650 ilast = ilast + 1 sofar = sofar + likfunc(ilast) if(sofar/intgrl.lt.1.d0-alpha.and.ilast.lt.nlikfunc)goto 650 if(sofar/intgrl.ge.1.d0-alpha.and. & ilast-ifirst.lt.ilbest-ifbest) then ifbest = ifirst ilbest = ilast endif if(ilast.lt.nlikfunc) goto 550 llimit = 1.d0+(ifbest+0.d0)/(unitdiv+0.d0) rlimit = 1.d0+(ilbest+0.d0)/(unitdiv+0.d0) return end c c EST1 c Routine which steps through the data sequence interval by interval, computing c the best tau, theta0, and theta1 for each. It returns taubest, t0best, and c t1best, the overall "winners" (overall maximum likelihood estimators). c subroutine est1(T,x,theta0,theta1,tau,lik,taubest,t0best, & t1best) double precision theta0(T), theta1(T), tau(T), lik(T) integer T, x(T) double precision loglbest,t1best,t0best,taubest,llik,sum0,sum1 double precision temp0, temp1 integer i, j double precision lnlik lik(1) = 0.d0 lik(2) = 0.d0 lik(T) = 0.d0 lik(T-1) = 0.d0 c Step through the data set interval by interval do 1000 i = 2,T-2 c Compute estimate for theta0(i) sum0 = 0.d0 do 100 j = 1, i sum0 = sum0 + x(j) 100 continue theta0(i) = sum0/(i+0.d0) c Compute estimate for theta1(i) sum1 = 0.d0 do 200 j = i+2, T sum1 = sum1 + x(j) 200 continue theta1(i) = sum1/(T-i-1.d0) c Compute estimate for tau(i) tau(i) = (x(i+1) + theta0(i)*i - theta1(i)*(i+1))/ & (theta0(i)-theta1(i)) c If the estimated change-point is not in the specified interval, c the MLE for the interval is one of the endpoints. if(tau(i).lt.i.or.tau(i).gt.i+1) then temp0 = 1.d0*(sum0 + x(i+1))/(i+1.d0) temp1 = 1.d0*(sum1 + x(i+1))/(T-i-0.d0) if(lnlik(T,i+0.d0,x,theta0(i),temp1).gt.lnlik(T, & i+1.d0,x,temp0,theta1(i))) then tau(i) = i theta1(i) = temp1 else tau(i) = i+1.d0 theta0(i) = temp0 endif endif llik=lnlik(T,tau(i),x,theta0(i),theta1(i)) c Compute log-likelihood for each i if(i.eq.2) then taubest = tau(i) t0best = theta0(i) t1best = theta1(i) loglbest = lnlik(T,tau(i),x,theta0(i),theta1(i)) else if(llik.gt.loglbest) then taubest = tau(i) t0best = theta0(i) t1best = theta1(i) loglbest = llik endif endif lik(i+1) = dexp(llik) 1000 continue return end c c ESTTHETA c Routine to compute estimates of theta0 and theta1 given a value for the c change-point tau. c subroutine esttheta(T,tau,x,theta0,theta1) double precision tau, theta0, theta1 integer T, x(T),sum double precision a, b, c if(a(T,tau).eq.0.d0) then theta1 = -c(T,tau,x)/b(T,tau,x) else theta1 = (-b(T,tau,x) + dsqrt(b(T,tau,x)**2 - 4.d0*a(T,tau)* & c(T,tau,x)))/(2.d0*a(T,tau)) endif sum = 0 do 100 i = 1,T sum = sum + x(i) 100 continue theta0 = (sum - (T-tau)*theta1)/tau if ( theta0 .le. 1.0d-13 .or. theta1 .le. 1.0d-13) then theta1 = (-b(T,tau,x) - dsqrt(b(T,tau,x)**2 - 4.d0*a(T,tau)* & c(T,tau,x)))/(2.d0*a(T,tau)) sum = 0 do 200 i = 1,T sum = sum + x(i) 200 continue theta0 = (sum - (T-tau)*theta1)/tau endif return end c c LNLIK c Routine to calculate the log-likelihood function for specified values of c tau, theta0, and theta1. c double precision function lnlik(T,tau,x,theta0,theta1) double precision tau, theta0, theta1, asum integer T,x(T),brack double precision p, zzlog lnlik = -theta0*brack(tau) asum = 0d0 do 100 i = 1,brack(tau) asum = asum + x(i) 100 continue lnlik = lnlik + zzlog(theta0)*asum asum = 0 do 200 i = brack(tau)+2,T asum = asum + x(i) 200 continue lnlik = lnlik + asum*zzlog(theta1) - theta1*(T-brack(tau)-1) lnlik = lnlik - p(tau)*theta0 - (1-p(tau))*theta1 & + zzlog(p(tau)*theta0 + (1-p(tau))*theta1)*x(brack(tau)+1) return end c c P c Function to give the fractional part of the number tau. c double precision function p(tau) double precision tau integer brack p = tau - brack(tau) return end c c c ZZLOG c Function to return the log of a number z, or a very large negative number if c the z is equal to zero. This is to avoid NaN's. c double precision function zzlog(z) double precision z if(z.eq.0.d0) then zzlog = -9.99d-100 else zzlog = dlog(z) endif return end c