



      real*8 pnew(2),aux,auxi,uni(4000000,4,4)
      real*8 err,thumb,dbetin
      real*8 xguess(12),xres(12),rp(5),fv(12),xsc(12),fsc(12)
      real*8 test,xnew(12),error,moy,std,test
      real*8 mu,sigma,xinf,xsup,aux1,aux2,moy1,moy2,std1,std2
      real*8 dnordf,dnorin
      integer ip(6),m,i,j,l,jj,t
      character rep

      real*8 erre,h,u(5000000,4,4)
      integer mc,nj,no,np

      external erset,thumb,drnun,drnnor,dbetin,dneqbf,rnset
      external dnordf,dnorin
      intrinsic dlog,dexp
      open(unit=60,file='data_multi_10_5.dat',status='new')

c data
      call rnset(642867)

      mc=100000
      nj=4
      no=2

c parameters of the beta distributions

      pnew(1)=3.0d0
      pnew(2)=3.0d0

c parameters of the truncated normal

       mu=3.0d0/4.0d0
       sigma=0.15d0

      xinf=0.0d0
      xsup=1.0d0
      aux1=dnordf((xinf-mu)/sigma)
      aux2=dnordf((xsup-mu)/sigma)-aux1
c      write(6,*)"aux1,aux2",aux1,aux2
c      read *,rep

c generate pseudo uniform random numbers

      do 47 m=1,mc
      do 47 i=1,nj
      do 47 j=1,no
        call drnun(1,uni(m,i,j))
   47 continue



      moy1=0.0d0
      std1=0.0d0
      moy2=0.0d0
      std2=0.0d0

      do 310 m=1,mc
        write(6,*)m
      do 310 i=1,nj
            aux=dbetin(uni(m,i,1),pnew(1),pnew(2))
            moy1=moy1+aux
            u(m,i,1)=aux

            aux=dnorin(uni(m,i,2)*aux2+aux1)
            u(m,i,2)=aux*sigma+mu
            moy2=moy2+u(m,i,2)
            if (u(m,i,2).lt.xinf.or.u(m,i,2).gt.xsup) then
               write(6,*)"problem with truncated normal",m,i,u(m,i,2)
               read *,rep
            endif
  310 continue


      moy1=moy1/(m*nj*1.0d0)
      moy2=moy2/(m*nj*1.0d0)
      do 315 m=1,mc
      do 315 i=1,nj
         std1=std1+(u(m,i,1)-moy1)**2
         std2=std2+(u(m,i,2)-moy2)**2
         write(60,100)u(m,i,1)
         write(60,100)u(m,i,2)
  315 continue
      std1=std1/(m*nj*1.0d0)
      std2=std2/(m*nj*1.0d0)

      write(6,*)"moy et std de t1",moy1,std1**0.5
      write(6,*)"moy et std de t2",moy2,std2**0.5

  100 format (2d20.12)

      end






