*...:....1....:....2....:....3....:....4....:....5....:....6....:....7....:....8
*    
*     mcarlo.f
*
*     Copyright (C) 1991.
*
*     Marie Davidian and A. Ronald Gallant 
*     P.O. Box 5513, Raleigh NC 27650-5513, USA   
*
*     Permission to use, copy, modify, and distribute this
*     software and its documentation for any purpose and
*     without fee is hereby granted, provided that the
*     above copyright notice appear in all copies and that
*     both that copyright notice and this permission notice
*     appear in supporting documentation.
*
*     This software is provided "as is" without any
*     expressed or implied warranty.
*
*...:....1....:....2....:....3....:....4....:....5....:....6....:....7....:....8
      subroutine mcarlo(mode,lvec,vec,objf,objgrd,nstate)
      implicit  real*8 (a-h,o-z)
      implicit  integer*4 (i-n)
      save

      include "global.f"

      real*8 objf,tmp
      real*8 vec,objgrd
      real*8 V,W,work
      real*8 sum,sumdel,u
      real*8 ftrial,fobs

      integer*4 mode,lvec,nstate
      integer*4 jseed,i,j,ii,jj,isw

      real*8 eps
      integer*4 mV,mwork

      parameter (mV=200,mwork=500,eps=1.d-300)

      dimension vec(lvec),objgrd(lvec)
      dimension V(mV),W(mV*mV),work(mwork)
      dimension sumdel(mtau+mtheta),u(mM)

*     Copy parameters NPSOL gets to see to vsav.

      j=0
      do i=1,lvsav
        if (ivec(i).eq.0) then
          j=j+1
          vsav(i)=vec(j)
        end if
      end do 

*     Make tau and theta agree with vsav.

      do i=1,ltau
        tau(i)=vsav(i)
      end do

      do i=1,ltheta
        theta(i)=vsav(ltau+i)
      end do

*     Error checks.

      if (ntermz.gt.mV) then
        write(3,3001) ntermz
        stop
      end if

*     If mwork is too small, whx should catch it.

      call whx(V,W,x,
     &  theta(iaA),theta(iaB0),theta(iaB),theta(iaR),
     &  M,L,ntermz,ntermx,Kz,Kx,work,mwork)

      ftrial=dfloat(npt)
      fobs=dfloat(nobs)
      jseed=iseed

      objf=0.d0

*     Anlalytic derivatives are not supplied for the polynomial
*     parameters of the SNP.  NPSOL computes them numerically.

      j=0
      do i=1,lvsav
        if (ivec(i).eq.0) then
          j=j+1
          if ((i.le.ltau).or.(i.ge.ltau+iaR)) then
            objgrd(j)=0.d0
          end if
        end if
      end do

*     Compute the log likelihood.

      do id=1,ncases

*     Integrate the d density with respect to the SNP density and
*     and store it in sum. 
*     Compute derivatives with respect to the parameters tau and 
*     R and store them in sumdel.

        sum=0.d0

        do j=1,lvsav
          sumdel(j)=0.d0
        end do

        do i=1,npt
          call sampz(jseed,u,V,W,M,ntermz,ntermx)

          do jj=1,M
            z(jj)=0.d0
          end do

          do jj=1,M
          do ii=1,jj
            z(ii)=z(ii)+theta(ibR+jj*(jj-1)/2+ii)*u(jj)
          end do
          end do

          isw=2
          call density(isw)

          sum=sum+d

          do jj=1,ltau
            sumdel(jj)=sumdel(jj)+dwtau(jj)
          end do

          do jj=1,M
          do ii=1,jj
            sumdel(ltau+ibR+jj*(jj-1)/2+ii)
     &        = sumdel(ltau+ibR+jj*(jj-1)/2+ii)
     &          + dwz(ii)*u(jj)
          end do
          end do
        end do

        sum = sum/ftrial

        do j=1,lvsav
          sumdel(j)=sumdel(j)/ftrial
        end do

        sum=sum+eps

        objf=objf-dlog(sum)/fobs

        tmp=sum*fobs

        j=0
        do i=1,lvsav
          if (ivec(i).eq.0) then
            j=j+1
            if ((i.le.ltau).or.(i.ge.ltau+iaR)) then
              objgrd(j)=objgrd(j)-sumdel(i)/tmp
            end if
          end if
        end do

      end do

      return
3001  format(1x,'error, mcarlo, increase mV to at least',i5)
      end
