double precision function zpdf(z) c****************************************************************** c c Function to evaluate the standard normal pdf at the c double precision value z c c******************************************************************* c double precision z c zpdf=dexp(-z*z/2.0d0)/dsqrt(8.0d0*datan(1.0d0)) c return end double precision function tpdf(t,ndf) c****************************************************************** c c Function to evaluate the pdf of the t distribution c having ndf degrees of freedom at the double precision value t. c c******************************************************************* c double precision v,half,one,two,vhalf,c1,c2,t,gama data half,one,two/0.5d0,1.0d0,2.0d0/ c v=dble(float(ndf)) vhalf=v/two c1=gama(half+vhalf)-gama(half)-gama(vhalf)-half*dlog(v) c2=(v+one)/two c tpdf=dexp(c1-c2*dlog(one+t*t/v)) c return end double precision function chipdf(chi,ndf) c****************************************************************** c c Function to evaluate the pdf of the chi-square distribution c having ndf degrees of freedom at the double precision value chi. c c******************************************************************* c double precision one,two,vhalf,c,chi,gama data one,two/1.0d0,2.0d0/ c if(chi.lt.0.0d0) then chipdf=0.0d0 return endif if(chi.eq.0.0d0) then chipdf=0.0d0 if(ndf.eq.2) chipdf=0.5d0 return endif vhalf=dble(float(ndf))/two c=-vhalf*dlog(two)-gama(vhalf) c chipdf=dexp(c-(chi/two)+(vhalf-one)*dlog(chi)) c return end double precision function fpdf(f,ndf1,ndf2) c****************************************************************** c c Function to evaluate the pdf of the F distribution having c degrees of freedom ndf1 and ndf2 at the double precision value f. c c******************************************************************* c double precision one,two,v1half,v2half,vsum,vrat,c1,f,gama data one,two/1.0d0,2.0d0/ c if(f.lt.0.0d0) then fpdf=0.0d0 return endif c v1half=dble(float(ndf1))/two v2half=dble(float(ndf2))/two vrat=v1half/v2half vsum=v1half+v2half c1=v1half*dlog(vrat)+gama(vsum)-gama(v1half)-gama(v2half) c if(f.eq.0.0d0) then fpdf=0.0d0 if(ndf1.eq.2) fpdf=dexp(c1) return endif c fpdf=dexp(c1+(v1half-one)*dlog(f)-vsum*dlog(one+vrat*f)) c return end subroutine zcdf(x,cdf,pdf) c*************************************************************** c c Subroutine to evaluate the cdf and pdf of the standard c normal distribution at the double precision value x using A&S, c 26.2.17 c c***************************************************************** c double precision x,one,pdf,t,t2,cdf,zpdf data one/1.0d0/ c pdf=zpdf(x) t=one/(one+.2316419d0*dabs(x)) t2=t*t cdf=one-pdf*(.319381530*t-.356563782*t2+1.781477937*t*t2- 1 1.821255978*t2*t2+1.330274429*t*t2*t2) if(x.lt.0.0d0) cdf=one-cdf c return end subroutine tcdf(x,nd,cdf,pdf) c****************************************************************** c c Subroutine to evaluate the cdf and pdf for the Student-t c distribution with nd degrees of freedom at the single point x c using the relationship of the t distribution with the c incomplete beta function. c c x (input) and cdf and pdf (output) are double precision c while nd is integer c c****************************************************************** c double precision zero,half,one,two,x,xx,bi,d,cdf,pdf,tpdf data zero,half,one,two/0.0d0,0.5d0,1.0d0,2.0d0/ c d=dble(float(nd)) xx=d/(d+x*x) call betai(xx,d/two,half,bi,ier) cdf=one-(bi/two) if(x.lt.zero) cdf=one-cdf pdf=tpdf(x,nd) c return end subroutine chicdf(x,nd,cdf,pdf) c******************************************************************* c c Subroutine to find the chi-square cdf and pdf for degrees of c freedom nd at the point x using the relationship of the c chi-square distribution with the incomplete gamma function. c c x (input) and cdf and pdf (output) are double precision c while nd (input) is integer c c c******************************************************************* c double precision zero,two,x,cdf,pdf,dhalf,xhalf,chipdf data zero,two/0.0d0,2.0d0/ c dhalf=dble(float(nd))/two xhalf=x/two call gammi(xhalf,dhalf,cdf,ier) pdf=chipdf(x,nd) c return end subroutine fcdf(x,nd1,nd2,cdf,pdf) c******************************************************************* c c Subroutine to find the F cdf and pdf for degrees of c freedom nd1 and nd2 at the point x using the relationship of c the F distribution with the incomplete beta function. c c x (input) and cdf and pdf (output) are double precision c while nd1 and nd2 (input) are integer c c c******************************************************************* c double precision two,x,cdf,pdf,v1,v2,xx,fpdf data two/2.0d0/ c v1=dble(float(nd1)) v2=dble(float(nd2)) xx=v2/(v2+v1*x) call betai(xx,v2/two,v1/two,cdf,ier) cdf=1.0d0-cdf if(x.le.0.0d0) cdf=0.0d0 pdf=fpdf(x,nd1,nd2) c return end subroutine zqnt(u,q,ier) c******************************************************************* c c Subroutine to evaluate the quantile function of the standard c normal distribution at the value u using A&S (26.2.22) as c starting value and Newton's method to find final value. c c u (input) and q (output) are double precision c ier (output) is an integer (0 means convergence, 1 means none) c c The calling routine must verify that u is between 0 and 1!! c c******************************************************************* c double precision u,u1,t,cdf,pdf,eps,del,q data itmax,eps/100,1.0d-10/ c ier=0 c c starting value: c u1=u if(u.gt.0.5d0) u1=1.0d0-u t=dsqrt(-2*dlog(u1)) q=(2.30753d0+.27061d0*t)/(1.0d0+.99229d0*t+.04481d0*t*t)-t if(u.gt.0.5d0) q=-q c c Iterate: c do 10 i=1,itmax call zcdf(q,cdf,pdf) if(pdf.lt.1.d-10) pdf=1.d-10 del=(cdf-u)/pdf if(dabs(cdf-u).lt.eps) go to 20 q=q-del 10 continue ier=1 20 continue return end subroutine tqnt(u,n,q,ier) c******************************************************************** c c Subroutine to evaluate the quantile function of the c t distribution with n degrees of freedom at the value u. c c u (input) and q (output) are double precision while n c (input) and ier (output) are integers c c ier = 0 (1) means convergence was (was not) achieved c c******************************************************************** c double precision eps,u,p,z1,z2,z3,z5,z7,z9,v1,v2,v3,v4, 1 q,pdf,cdf,del c itmax=100 eps=1.0d-10 ier=0 c c Find Starting Value as in A \& S, 26.7.5: c p=u if(u.le.0.5d0) p=1.0d0-u call zqnt(p,z1,ier) z2=z1*z1 z3=z1*z2 z5=z2*z3 z7=z5*z2 z9=z7*z2 q=z1 v1=dble(float(n)) v2=v1*v1 v3=v1*v2 v4=v2*v2 q=q+(z3+z1)/(4.0d0*v1) q=q+(5.0d0*z5+16.0d0*z3+3.0d0*z1)/(96.0d0*v2) q=q+(3.0d0*z7+19.0d0*z5+17.0d0*z3-15.0d0*z1)/(384.0d0*v3) q=q+(79.0d0*z9+776.0d0*z7+1482.0d0*z5-1920.0d0*z3- 1 945.0d0*z)/(92160.0d0*v4) c c Iterate: c do 10 i=1,itmax call tcdf(q,n,cdf,pdf) if(pdf.lt.1.d-10) pdf=1.d-10 del=(cdf-p)/pdf if(dabs(cdf-p).lt.eps) go to 20 q=q-del 10 continue ier=1 20 continue if(u.le.0.5d0) q=-1.0d0*q 99 return end subroutine chiqnt(u,n,q,ier) c******************************************************************** c c Subroutine to evaluate the quantile function of the c chi-square distribution with n degrees of freedom at the value u. c c u (input) and q (output) are double precision while n c (input) and ier (output) are integers c c ier = 0 (1) means convergence was (was not) achieved c c******************************************************************** double precision two,eps,u,y,v,chk,q,cdf,pdf,del,gama c c itmax=100 eps=1.0d-10 ier=0 c c Find Starting Value as on pg 118 of Kennedy and Gentle: c call zqnt(u,y,ier) two=2.0d0 v=dble(float(n)) chk=-1.24d0*dlog(u) if(v.lt.chk) then q=dexp((two/v)*(gama(v/two)+dlog(u*v)+ 1 ((v-two)/two)*dlog(two))) else q=dexp(dlog(v)+3.0d0*dlog(1.0d0-(two/(9.0d0*v))+ 1 y*dexp(0.5d0*dlog(two/(9.0d0*v))))) if(q.gt.2.20d0*v+6.0d0) then q=-two*(dlog(1.0d0-u)- 1 ((v/two)-1.0d0)*dlog(q/two)+gama(v/two)) endif endif c c Iterate: c do 10 i=1,itmax call chicdf(q,n,cdf,pdf) if(pdf.lt.1.d-10) pdf=1.d-10 del=(cdf-u)/pdf if(dabs(cdf-u).lt.eps) go to 20 q=q-del 10 continue ier=1 20 continue 99 continue return end subroutine fqnt(u,n,m,q,ier) c******************************************************************** c c Subroutine to evaluate the quantile function of the c F(n,m) distribution at the value u. c c u (input) and q (output) are double precision while n and m c (input) and ier (output) are integers c c ier = 0 (1) means convergence was (was not) achieved c c******************************************************************** double precision a,b,d,e,h,p,q,t,eps,u,w,x,x1,x2,y,y1,y2 double precision pdf,cdf,del c c itmax=100 eps=1.0d-10 ier=0 c c handle n=1 or m=1 case: c if(n.eq.1) then h=(1.0d0+u)/2.0d0 call tqnt(h,m,q,ier) q=q**2 go to 99 endif if(m.eq.1) then h=1.d0-(u/2.0d0) call tqnt(h,n,q,ier) q=1.0d0/(q*q) go to 99 endif c c handle case where both n and m are greater than 1: c c Find Starting Value: c call zqnt(u,y,ier) c c 26.5.22 and 26.6.16: c a=2.0d0*dble(float(m))-1.0d0 b=2.0d0*dble(float(n))-1.0d0 h=2.0d0/((1.0d0/a)+(1.0d0/b)) d=(y*y-3.0d0)/6.0d0 w=(y*dsqrt(h+d)/h)- 1 ((1.0d0/b)-(1.0d0/a))*(d+(5.0d0/6.0d0)-(2.0d0/(3.0d0*h))) q=dexp(2.0d0*w) c c Iterate: c do 10 i=1,itmax call fcdf(q,n,m,cdf,pdf) if(pdf.lt.1.d-10) pdf=1.d-10 del=(cdf-u)/pdf if(dabs(cdf-u).lt.eps) go to 20 q=q-del 10 continue ier=1 20 continue 99 continue return end double precision function gama(x) c***************************************************************** c c Function to find the log of the gamma function as on page 157 c of Numerical Recipes. c c x and gama are double precision c c calling routine must verify that x is positive!!! c c******************************************************************* c double precision cof(6),stp,half,one,fpf,x,tmp,ser,xx data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0, 1 -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/ data half,one,fpf/0.5d0,1.0d0,5.5d0/ c xx=x-one tmp=xx+fpf tmp=(xx+half)*dlog(tmp)-tmp ser=one do 10 j=1,6 xx=xx+one 10 ser=ser+cof(j)/xx gama=tmp+dlog(stp*ser) c return end subroutine gammi(x,a,gi,ier) c***************************************************************** c c Subroutine to find the incomplete gamma function as on page 161 c of Numerical Recipes. c c x and a (input) and gi (output) are double precision c ier (output) is 0 if no error, and 1 if an error c c calling routine must verify that x and a are positive!! c c******************************************************************* c real*8 x,a,gi,one,gama,gln,ap,sum,del,eps,zero,gold, 1 a0,b0,a1,b1,fac,an,ana,anf,g data zero,one,eps,itmax/0.0d0,1.0d0,1.0d-8,100/ c ier=0 if(x.le.zero) then gi=zero go to 99 endif gln=gama(a) c c handle x.lt.a+1: Infinite Series c if(x.lt.a+one) then ap=a sum=one/a del=sum do 10 n=1,itmax ap=ap+one del=del*(x/ap) sum=sum+del if(dabs(del).lt.dabs(sum)*eps) then gi=sum*dexp(-x+a*dlog(x)-gln) go to 99 endif 10 continue ier=1 go to 99 else c c Handle x.ge.a+1: Continued Fraction c gold=zero a0=one a1=x b0=zero b1=one fac=one do 20 n=1,itmax an=float(n) ana=an-a a0=(a1+a0*ana)*fac b0=(b1+b0*ana)*fac anf=an*fac a1=x*a0+anf*a1 b1=x*b0+anf*b1 if(a1.ne.zero) then fac=one/a1 g=b1*fac if(dabs((g-gold)/g).lt.eps) then gi=one-dexp(-x+a*dlog(x)-gln)*g go to 99 endif gold=g endif 20 continue ier=1 go to 99 c endif c 99 continue return end subroutine betai(x,a,b,bi,ier) c***************************************************************** c c Subroutine to find the incomplete beta function as on page 167 c of Numerical Recipes. c c x, a, and b (input) and bi (output) are double precision c ier (output) is 0 if no error, and 1 if an error c c Calling routine must verify that x is between 0 and 1 and that c a and b are positive!! c c******************************************************************* c double precision x,a,b,gama,bi,eps,one,two,bt,x1,a1,b1,am,bm, 1 az,qab,qap,qam,bz,em,tem,d,ap,bp,app,bpp,aold data one,two,eps,itmax/1.0d0,2.0d0,1.0d-10,100/ c ier=0 if(x.le.0.0d0) then bi=0.0d0 go to 99 endif if(x.ge.one) then bi=one go to 99 endif bt=dexp(gama(a+b)-gama(a)-gama(b)+a*dlog(x)+b*dlog(one-x)) x1=x a1=a b1=b isym=0 if(x.ge.(a+one)/(a+b+two)) then isym=1 a1=b b1=a x1=one-x endif c am=one bm=one az=one qab=a1+b1 qap=a1+one qam=a1-one bz=one-qab*x1/qap do 10 m=1,itmax em=m tem=em+em d=em*(b1-em)*x1/((qam+tem)*(a1+tem)) ap=az+d*am bp=bz+d*bm d=-(a1+em)*(qab+em)*x1/((a1+tem)*(qap+tem)) app=ap+d*az bpp=bp+d*bz aold=az am=ap/bpp bm=bp/bpp az=app/bpp bz=one if(dabs(az-aold).lt.eps*dabs(az)) go to 20 10 continue ier=1 go to 99 20 if(isym.eq.1) bi=one-bt*az/b if(isym.eq.0) bi=bt*az/a c 99 continue return end subroutine runif(dseed,n,x) c**************************************************************** c c Subroutine to fill the first n elements of the double precision c array x with iid U(0,1) random variables. c c The argument dseed is double precision and on input contains c a large whole number for a seed. On output it contains a new c seed that can be used as input the next time runif is called. c c**************************************************************** c double precision x(n) real*8 dseed,svnto5,t31m1 data svnto5,t31m1/16807.d0,2147483647.d0/ c do 10 i=1,n dseed=dmod(svnto5*dseed,t31m1) 10 x(i)=dseed/t31m1 c return end subroutine rnorm(dseed,n,x) c***************************************************************** c c Subroutine to fill the first n elements of the double precision c array x with iid N(0,1) random numbers. c c The argument dseed is double precision and on input contains c a large whole number for a seed. On output it contains a new c seed that can be used as input the next time rnorm is called. c c******************************************************************* c double precision x(n) double precision dseed,v,y,u(2) i=0 5 call runif(dseed,2,u) u(1)=2*(u(1)-.5) u(2)=2*(u(2)-.5) v=u(1)*u(1)+u(2)*u(2) if(v.ge.1.0) go to 5 y=sqrt(-2.*log(v)/v) i=i+1 x(i)=u(1)*y if(i.eq.n) return i=i+1 x(i)=u(2)*y if(i.eq.n) return go to 5 return end subroutine qsort(x,n) c********************************************************************* c c Quicksort a real array x of length n. c c******************************************************************** c double precision x(n) integer ll(20),lr(20) c c ns is # of pieces left to partition c c ll and lr are left and right borders of pieces c ns=1 ll(1)=1 lr(1)=n c c After splitting we come here. Stop when no pieces left. c 5 if(ns.eq.0) go to 99 c i=ll(ns) j=lr(ns) nl=j-i+1 c c Sort pieces of size 1 or 2: c if(nl.le.2) then ns=ns-1 if(nl.eq.2.and.x(i).gt.x(j)) call swap(x(i),x(j)) go to 5 endif c c Sort first, ``middle'' and last elements: c nm=(i+j)/2 if(x(i).gt.x(nm)) call swap(x(i),x(nm)) if(x(nm).gt.x(j)) call swap(x(nm),x(j)) if(x(i).gt.x(nm)) call swap(x(i),x(nm)) c c If piece of size 3, it's now sorted c if(nl.eq.3) then ns=ns-1 go to 5 endif c c Put middle (target) into 1st element and keep a copy for comparisons c ax=x(nm) call swap(x(i),x(nm)) c c Now we'll look for 1st one from left (starting with 2nd) > target c and first from right < target. If no such pair we end up at 20. c itemp=i i=i+1 c 10 if(x(i).le.ax) then if(i.eq.j) go to 20 i=i+1 go to 10 endif c 15 if(x(j).ge.ax) then if(i.eq.j) go to 20 j=j-1 go to 15 endif c c Can only get to here if we found a pair to swap c call swap(x(i),x(j)) c c Only go back to look for another pair if there's a chance to get one c if(j-i.le.1) then go to 20 else i=i+1 j=j-1 go to 10 endif c c Get to here when no more pairs to swap. i might be 1 space too far c to right c 20 if(x(i).gt.ax) i=i-1 call swap(x(itemp),x(i)) c c Now put shorter piece at end of ll and lr (this guarantees small c number of pieces at any one time). A piece might be empty, but that c is taken care of above. c nleft=i-ll(ns) nright=lr(ns)-i ns=ns+1 if(nleft.le.nright) then ll(ns)=ll(ns-1) ll(ns-1)=i+1 lr(ns)=i-1 else lr(ns)=lr(ns-1) lr(ns-1)=i-1 ll(ns)=i+1 endif go to 5 c 99 continue return end subroutine swap(a,b) c********************************************************************* c c Swap two reals. c c******************************************************************** c double precision a,b,c c=b b=a a=c c return end subroutine MGS(X,n,m,ndim,mdim,R,d,ier) c******************************************************************* c c Subroutine to find Modified Gram-Scmidt Decomposition X=QR of an c n by m matrix X. c c Input: X, n, m c ndim, mdim: row dimension of X and R in calling program c c Output: X: contains Q c R c d: vector containing diagonal elements of Q transpose Q c ier: 0 (1) means X is (is not) nonsingular c c Subprograms called: double precision function dbprod c c******************************************************************* double precision X(ndim,1),R(mdim,1),d(1) double precision dbprod,eps data eps/1.e-20/ ier=1 do 100 i=1,m r(i,i)=1.0d0 d(i)=dbprod(X(1,i),X(1,i),n) if(d(i).lt.eps) go to 99 if(i.eq.m) then ier=0 go to 99 endif do 30 j=i+1,m alpha=dbprod(X(1,j),X(1,i),n)/d(i) R(i,j)=alpha R(j,i)=0.0d0 do 20 k=1,n 20 X(k,j)=X(k,j)-alpha*X(k,i) 30 continue 100 continue 99 continue return end double precision function dbprod(x,y,n) c---------------------------------------------------------------- c c c------------------------------------------------------------------- double precision x(n),y(n) dbprod=0.0d0 do 10 i=1,n 10 dbprod=dbprod+x(i)*y(i) return end subroutine solvch(a,b,n,ndim,al,d,wk,x,ier) c********************************************************************** c c Subroutine to solve the system of equations a*x=b, where a is an c n by n symmetric matrix that is supposed to be positive definite. c c ndim is the row dimension of a and the n by n matrix al in the c calling routine. c c d is a work vector of length n c c ier is an output error indicator; 0 means no error, anything else means c a is not positive definite. c c subprograms called: mchol c c*********************************************************************** c double precision a(ndim,1),b(n),x(n),al(ndim,1),d(n),wk(n) double precision c c call mchol(a,n,ndim,al,d,ier) if(ier.ne.0) return do 20 i=1,n c=b(i) if(i.gt.1) then do 10 j=1,i-1 10 c=c-al(i,j)*wk(j) endif 20 wk(i)=c do 30 i=1,n 30 wk(i)=wk(i)/d(i) do 50 i=n,1,-1 c=wk(i) if(i.lt.n) then do 40 j=i+1,n 40 c=c-al(i,j)*x(j) endif 50 x(i)=c end subroutine mchol(a,n,ndim,al,d,ier) c************************************************************************ c c Subroutine to obtain the Cholesky factors al (L in lower triangle, c L' in upper triangle) and d (diagonal elements of D in the c vector d) of the (n by n) positive definite matrix a, ie a=L*D*L'. c c ndim is the row dimension of a and al in the calling routine. c c ier is an output error indicator; 0 means no error, anything else c means a is not positive definite. c c subprograms called: none c c************************************************************************ c double precision a(ndim,1),al(ndim,1),d(1) double precision c data eps/1.e-20/ c Initialize: ier=0 al(1,1)=1.0 d(1)=a(1,1) if(d(1).lt.eps) then ier=1 return endif if(n.eq.1) return do 10 i=2,n c Get ith row of L: al(i,i)=1.0 do 20 j=1,i-1 c=dble(a(i,j)) if(j.gt.1) then do 30 l=1,j-1 30 c=c-dble(al(i,l))*dble(d(l))*dble(al(j,l)) endif al(i,j)=c/d(j) 20 al(j,i)=al(i,j) c Get D(i,i): c=dble(a(i,i)) do 50 l=1,i-1 50 c=c-dble(d(l))*dble(al(i,l))**2 d(i)=c c Check for positivity of D(i,i): if(d(i).lt.eps) then ier=i return endif 10 continue return end subroutine sweep(a,mdim,m,k1,k2,ier) c --------------------------------------------------------------- c c Subroutine to sweep the mxm matrix a on its k1 st thru k2 th c diagonals. ier is 1 if a is singular. mdim is row dimension c of a in calling routine. c c --------------------------------------------------------------- double precision a(mdim,1) double precision d ier=1 do 50 k=k1,k2 if(abs(a(k,k)).lt.1.e-20) return d=1./a(k,k) a(k,k)=1. do 10 i=1,m a(k,i)=d*a(k,i) 10 if(i.ne.k) a(i,k)=-a(i,k)*d do 20 i=1,m do 20 j=1,m 20 if((i.ne.k).and.(j.ne.k)) a(i,j)=a(i,j)+a(i,k)*a(k,j)/d 50 continue ier=0 return end subroutine nr(theta,n,ndim,maxit,eps,grdhss,g,h,al,d,wk, 1 del,ier) c--------------------------------------------------------------- c c Subroutine to perform Newton-Raphson minimization of a c function of n variables. c c Input: c theta: vector of length n with starting values c n c ndim: row dimension of h and al in calling routine c maxit: maximum number of iterations to perform c eps: convergence criterion, how close successive c iterations must be to conclude that c convergence has been reached. c grdhss: name of the subroutine that will calculate the c gradient and Hessian (must be declared c external in calling routine). c c Work: c al: (n by n) matrix having ndim rows c d, wk, del: vectors of length n c c Output: c theta: value of variables at last iteration c g, h: gradient and Hessian at final iteration c ier: error indicator (0 means none, -1 means no c convergence, + means H not positive definite). c c Subprograms Called: grdhss, solvch (mchol) c c--------------------------------------------------------------- double precision theta(n),g(n),h(ndim,1),al(ndim,1),d(n), 1 wk(n),del(n) double precision eps do 100 it=1,maxit call grdhss(theta,n,ndim,g,h) call solvch(h,g,n,ndim,al,d,wk,del,ier) if(ier.ne.0) return do 10 i=1,n 10 if(abs(del(i)).ge.(eps*abs(theta(i)))) go to 20 ier=0 return 20 if(it.ge.maxit) then ier=-1 return endif do 30 i=1,n 30 theta(i)=theta(i)-del(i) 100 continue return end