        subroutine ghr(x,y,sig2,b,iseed)
        implicit real*8 (a-h,o-z)
        real*8 x(10,250), y(250), b(10)
        real*8 sig2
        integer iseed
        call Gbeta(x,y,sig2,b,iseed)
        call Gsig2(x,y,b,sig2,iseed)
        return
        end



        subroutine Gbeta(x,y,sig2,b,iseed)
        implicit real*8 (a-h,o-z)
        real*8 x(10,250), y(250), b(10)
        real*8 sig2
        real*8 sum1,sum2,amu
        integer iseed
        real*8 aleft,aright,amean,sigma,tn,psi
        logical la,lb
        la=.false.
        lb=.false.
c      
        do 100 k=1,10
         sum1=0.0d0
         sum2=0.0d0
         do 110 i=1,207
          sum1=sum1+x(k,i)**2
          amu=0.0d0
          do 120 j=1,10
           if (j .eq. k) goto 120
           amu=amu+b(j)*x(j,i) 
 120      continue
          sum2=sum2+(y(i)-amu)*x(k,i)
 110     continue
         amean=sum2/sum1
         sigma=dsqrt( sig2/sum1 )
         if (k .eq. 1) then
          aright=(b(2)-amean)/sigma
          aleft=(-amean)/sigma
         endif
         if ( (k .gt. 1) .and. (k .lt. 10) ) then
          aleft=(b(k-1)-amean)/sigma
          aright=(b(k+1)-amean)/sigma
         endif
         if ( k .eq. 10 ) then
          aleft=(b(9)-amean)/sigma
          lb=.true.
          psi=(.21110**2*2.0d0)/(.21110**2*2.0d0+
     1         sigma**2)
          amean=psi*amean+(1.0d0-psi)*0.99797
          sigma=dsqrt( (1.0d0-psi)*(.21110**2*2.0d0) )
         endif
         tn=YTUVN(aleft,aright,la,lb,iseed)
         lb=.false.
         la=.false.
         b(k)=amean+tn*sigma
 100    continue
c
        return
        end


       subroutine Gsig2(x,y,b,sig2,iseed)
       integer n,iseed
       implicit real*8 (a-h,o-z)
       real*8 x(10,250),y(250),b(10)
       real*8 sig2,amu,alam(250)
       real*8 scaler,shape,r(1)
       external  DRNGAM
       shape=207.0d0/2.0d0
       scaler=0.0d0
       do 220 i=1,207
        amu=0.0d0
        do 230 j=1,10
         amu=amu+b(j)*x(j,i)
 230    continue
        scaler=scaler+(y(i)-amu)**2
 220   continue
       scaler=scaler/2.0d0
       call rnset( iseed )
       call DRNGAM(1,shape,r)
       call rnget( iseed )
       sig2=scaler/r(1)
       return
       end
