PROGRAM spect_mix_logit_rj_sim2
USE lin_eig_self_int
USE show_int
USE wrrrn_int
USE rnset_int
IMPLICIT none
INTERFACE
FUNCTION log_spect_ar2(nfreq,comp)
IMPLICIT none
INTEGER, INTENT(IN) :: nfreq, comp
DOUBLE PRECISION, DIMENSION(nfreq) :: log_spect_ar2
END FUNCTION
END INTERFACE
INTEGER, PARAMETER :: nobs=64, nbasis=10
INTEGER, PARAMETER :: nchunks=16
INTEGER, PARAMETER :: nsim=100
INTEGER, PARAMETER :: nfreq=nobs/2 + 1
INTEGER, PARAMETER :: ny=nchunks*nfreq
INTEGER, PARAMETER :: nbeta=nbasis+1
INTEGER, PARAMETER :: nloop=100000,nwarmup=20000
INTEGER :: i, j, k,iseed
DOUBLE PRECISION, allocatable :: tp(:),omega(:,:),x(:,:)
DOUBLE PRECISION, allocatable :: xx(:,:),y(:,:),yy(:,:),xdelta(:,:)
DOUBLE PRECISION, allocatable :: d(:), v_s(:,:)
DOUBLE PRECISION :: sigalpha,sigdelta
DOUBLE PRECISION, DIMENSION(nfreq,nchunks) :: fit
DOUBLE PRECISION, DIMENSION(nfreq,nchunks) :: log_spect_true
INTEGER, DIMENSION(nchunks,nsim) :: true_ind
DOUBLE PRECISION :: mse
DOUBLE PRECISION, DIMENSION(nchunks-1) :: clevel
INTEGER, DIMENSION(nchunks-1) :: iclson,icrson
open(10, file = 'log_per_ar10_100sim')
open(12, file = 'iter_ar10_100sim22')
open(15, file = 'fitted_ar10_100sim22')
open(20, file = 'lcl_ar10_100sim22')
open(22, file = 'ucl_ar10_100sim22')
iseed=24378
!iseed=54674
call rnset(iseed)
ALLOCATE(tp(nfreq))
ALLOCATE(omega(nfreq,nfreq))
DO i = 1, nfreq
!tp(i)=(dble(i)-1.0d0)/(2.0d0*(dble(nfreq)-1.0d0))
tp(i)=(2.0d0*dble(i)-1.0d0)/(2.0d0*dble(nfreq))
END DO
DO i=1,nfreq
DO j=i,nfreq
omega(i,j)=tp(i)**2*(tp(j)-tp(i)/3.0d0)/2.0d0
END DO
END DO
DO i=2,nfreq
DO j=1,i-1
omega(i,j)=omega(j, i)
END DO
END DO
!call wrrrn('omega',omega,nra=nfreq)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ALLOCATE(d(nfreq),v_s(nfreq,nfreq))
!Find the eigenvalue decomposition
call lin_eig_self(omega,d,v=v_s)
ALLOCATE(x(nfreq,nbasis))
DO i=1,nfreq
DO j=1,nbasis
x(i,j)=(-1)**(j+1)*v_s(i,j)*sqrt(d(j))
END DO
END DO
!Note: multiply by (-1)**(j+1) to match matlab's eigenvectors
!do i=1,nfreq
! print*, x(i,3)
!end do
!construct the matrix xx in which the first column is 1's
!and the rest of the columns are the matrix x. This combines
!the Z and X matrices. Z is just a column of 1's since alpha_1
!is constrained to be zero
ALLOCATE(xx(nfreq,nbeta))
xx(:,1)=1.0d0
DO j=2,nbeta
xx(:,j)=x(:,nbasis-j+2)
END DO
!DO j=1,nbeta
! DO i=1,nfreq
! read(8,*) xx(i,j)
! END DO
!END DO
DEALLOCATE(x)
ALLOCATE(xdelta(nchunks,2))
xdelta(:,1)=1.0d0
xdelta(:,2)=dble((/(k,k=1,nchunks)/))/dble(nchunks)
sigalpha=100.0d0
sigdelta=4.0d0
ALLOCATE(yy(nchunks*nfreq,nsim),y(nfreq,nchunks))
DO j=1,nsim
DO i=1,nchunks*nfreq
read(10,*) yy(i,j)
END DO
END DO
! DO i=1,nchunks
! IF(i.LE.8) true_ind(i,:)=1
! IF(i.GE.9 .AND. i.LE.12) true_ind(i,:)=2
! IF(i.GE.13 .AND. i.LE.16) true_ind(i,:)=3
! END DO
DO i=1,nchunks
IF(i.LE.3) true_ind(i,:)=1
IF(i.GE.4 .AND. i.LE.16) true_ind(i,:)=2
END DO
DO k=98,98
print*,k
y=reshape(yy(:,k),(/nfreq,nchunks/))
call revjump(iseed,nloop,nwarmup,nfreq,nchunks,nbeta,&
nbasis,sigalpha,sigdelta,xx,xdelta,y,fit,&
clevel,iclson,icrson)
!DO j=1,nchunks
! log_spect_true(:,j)=log_spect_ar2(nfreq,true_ind(j,k))
!END DO
!mse=sum((fit-log_spect_true)**2)/dble(nchunks*nfreq)
DO j=1,nchunks
DO i=1,nfreq
write(15,*) fit(i,j)
END DO
END DO
!write(20,*) mse
END DO
END PROGRAM spect_mix_logit_rj_sim2
!
!
SUBROUTINE revjump(iseed,nloop,nwarmup,nfreq,nchunks,nbeta,&
nbasis,sigalpha,sigdelta,xx,xdelta,y,fit,clevel,iclson,icrson)
!USE show_int
USE wrrrn_int
USE rnunf_int
USE linear_operators
USE lin_eig_self_int
USE eqtil_int
USE maxb_m, ONLY: maxb
IMPLICIT none
INTERFACE
FUNCTION findloc(m,n,ind,ii)
INTEGER, INTENT(IN) :: m,n
INTEGER, DIMENSION(n), INTENT(IN) :: ind
INTEGER, DIMENSION(m) :: findloc
INTEGER :: ii,j,k
END FUNCTION
FUNCTION mvrnorm(iseed,n,mean,var)
USE lin_eig_self_int
USE linear_operators
IMPLICIT none
INTEGER, INTENT(IN) :: n,iseed
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: mean
DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var
DOUBLE PRECISION, DIMENSION(n) :: mvrnorm
END FUNCTION
FUNCTION mvrnorm1(iseed,n,mean,var)
USE linear_operators
USE show_int
IMPLICIT none
INTEGER, INTENT(IN) :: n,iseed
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: mean
DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var
DOUBLE PRECISION, DIMENSION(n) :: mvrnorm1
END FUNCTION
FUNCTION log_like(nfreq,nchunks,nbeta,nexp,xx,y,beta,prior_ind)
IMPLICIT none
INTEGER, INTENT(IN) :: nfreq,nchunks,nbeta,nexp
DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: xx
DOUBLE PRECISION, DIMENSION(nfreq,nchunks), INTENT(IN) :: y
DOUBLE PRECISION, DIMENSION(nexp,nchunks), INTENT(IN) :: prior_ind
DOUBLE PRECISION, DIMENSION(nexp,nbeta), INTENT(IN) :: beta
DOUBLE PRECISION :: log_like
END FUNCTION
FUNCTION cumsum(n,x)
IMPLICIT none
INTEGER, INTENT(IN) :: n
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x
DOUBLE PRECISION, DIMENSION(n) :: cumsum
END FUNCTION
END INTERFACE
INTEGER, INTENT(IN) :: nfreq,nchunks,nbeta,nbasis,nloop,&
nwarmup,iseed
DOUBLE PRECISION, INTENT(IN) :: sigalpha,sigdelta
DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: xx
DOUBLE PRECISION, DIMENSION(nchunks,2), INTENT(IN) :: xdelta
DOUBLE PRECISION, DIMENSION(nfreq,nchunks), INTENT(IN) :: y
DOUBLE PRECISION, DIMENSION(nfreq,nchunks), INTENT(OUT) :: fit
INTEGER, PARAMETER :: nquan=2, nexp_max=10
INTEGER :: nexperts
DOUBLE PRECISION, ALLOCATABLE :: beta(:,:),delta(:,:),tau(:),prior_ind(:,:),&
prior_ind_prop(:,:),delta_prop(:,:),tau_prop(:),&
beta_prop(:,:)
DOUBLE PRECISION, DIMENSION(2) :: mean_d
DOUBLE PRECISION, DIMENSION(2,2) :: var_d
DOUBLE PRECISION, DIMENSION(nbeta) :: mean_b
DOUBLE PRECISION, DIMENSION(nbeta,nbeta) :: var_b
DOUBLE PRECISION, ALLOCATABLE :: h(:,:)
DOUBLE PRECISION :: u,drnunf
DOUBLE PRECISION :: u1,u2
DOUBLE PRECISION :: log_like_curr,log_like_new,met_rat_jump,epsilon_jump
INTEGER :: ndelta,nproposed,counter
INTEGER :: i,j,k,iter
INTEGER, DIMENSION(nchunks) :: ind0
INTEGER, DIMENSION(nchunks) :: ind,kk
INTEGER, ALLOCATABLE :: kk_ind(:)
INTEGER :: n_chunk_ind
DOUBLE PRECISION, ALLOCATABLE :: hh(:)
DOUBLE PRECISION, ALLOCATABLE :: post_ind(:,:)
DOUBLE PRECISION, DIMENSION(nbeta) :: param
DOUBLE PRECISION, DIMENSION(nfreq) :: fit_chunk
DOUBLE PRECISION, DIMENSION(nfreq,nchunks) :: fit_ch
DOUBLE PRECISION, DIMENSION(nchunks) :: sum_fit_ch
DOUBLE PRECISION, DIMENSION(nfreq,nchunks) :: fhat
DOUBLE PRECISION, DIMENSION(nfreq,nchunks) :: cum_fhat
DOUBLE PRECISION, DIMENSION(nchunks-1) :: dist
DOUBLE PRECISION, DIMENSION(nchunks-1) :: clevel
INTEGER, DIMENSION(nchunks-1) :: iclson,icrson
DOUBLE PRECISION, DIMENSION(nfreq,3) :: fit_exp,fit_exp_sum
INTEGER :: count_fit, nmiss
DOUBLE PRECISION :: qprop(nquan), quan(nquan), xlo(nquan), xhi(nquan)
DOUBLE PRECISION :: lcl(nfreq,nchunks), ucl(nfreq,nchunks)
DOUBLE PRECISION :: fit_it(nfreq,nloop-nwarmup,nchunks)
!DOUBLE PRECISION :: prior_ind_hat(4,nchunks)
INTEGER :: acc_death, acc_birth
!
nexperts=1
ndelta=2*(nexperts-1)
ALLOCATE(tau(nexperts),beta(nexperts,nbeta),delta(2,nexperts),&
prior_ind(nexperts,nchunks),h(nexperts,nchunks))
tau(1)=10000.0d0
!Note that if nexperts<2 the loop is not executed
DO j=2,nexperts
tau(j)=tau(j-1)/10.0d0
END DO
delta=0.0d0
mean_d=0.0d0
mean_b=0.0d0
var_d=sigdelta*eye(2)
!call em(iseed,nfreq,nchunks,nbeta,ndelta,nexperts,nbasis,sigalpha,&
! sigdelta,xx,xdelta,y,ind0,beta,prior_ind)
!ind=ind0
DO j=1,nexperts
call drnun(nchunks,h(j,:))
END DO
ind=maxloc(h,dim=1)
DEALLOCATE(h)
DO j=1,nexperts
kk=0
where(ind.eq.j) kk=1
n_chunk_ind=sum(kk)
ALLOCATE(kk_ind(n_chunk_ind),hh(n_chunk_ind))
!find the chunks belonging to expert j
kk_ind=findloc(n_chunk_ind,nchunks,ind,j)
hh=1.0d0
param=0.0d0
call maxb(param,tau(j),hh,sigalpha,nfreq,n_chunk_ind,nbeta,&
nexperts,nbasis,xx,y(:,kk_ind),beta(j,:))
DEALLOCATE(kk_ind,hh)
END DO
fit=0.0d0
!!!!!!!!!!!!!!!!
count_fit=0
fit_exp_sum=0.0d0
!!!!!!!!!!!!!!!
acc_birth=0
acc_death=0
DO iter=1,nloop
!print*,'iter',iter
!print*,'nexperts',nexperts
write(12,*) iter,nexperts
!dimension jumping step
u=rnunf()
IF(u>0.5d0 .OR. nexperts==1) THEN
!BIRTH STEP
nproposed=min(nexp_max,nexperts+1)
ALLOCATE(beta_prop(nproposed,nbeta),delta_prop(2,nproposed),&
prior_ind_prop(nproposed,nchunks),tau_prop(nproposed))
!generating new parameters beta, delta and tau
DO j=1,nexperts
delta_prop(:,j)=delta(:,j)
beta_prop(j,:)=beta(j,:)
tau_prop(j)=tau(j)
END DO
delta_prop(:,nproposed)=mvrnorm1(iseed,2,mean_d,var_d)
tau_prop(nproposed)=drnunf()*tau(nexperts)
var_b=tau_prop(nproposed)*eye(nbeta)
var_b(1,1)=sigalpha
beta_prop(nproposed,:)=mvrnorm1(iseed,nbeta,mean_b,var_b)
DO j=1,nproposed
prior_ind_prop(j,:)=exp(matmul(xdelta,delta_prop(:,j)))/&
sum(exp(matmul(xdelta,delta_prop)),2)
END DO
DO j=1,nexperts
prior_ind(j,:)=exp(matmul(xdelta,delta(:,j)))/&
sum(exp(matmul(xdelta,delta)),2)
END DO
!Checking dominance condition
counter=0
DO j=1,nproposed
DO k=1,nproposed
IF(maxval(prior_ind_prop(j,:)-prior_ind_prop(k,:))<0.0d0) &
counter=counter+1
END DO
END DO
log_like_curr=log_like(nfreq,nchunks,nbeta,nexperts,xx,y,beta,prior_ind)
log_like_new=log_like(nfreq,nchunks,nbeta,nproposed,xx,y,beta_prop,&
prior_ind_prop)
met_rat_jump=log_like_new-log_like_curr
IF(nexperts==1) THEN
epsilon_jump=min(1.0d0,0.5d0*exp(met_rat_jump))
ELSE IF (nexperts==nexp_max-1) THEN
epsilon_jump=min(1.0d0,2.0d0*exp(met_rat_jump))
ELSE
epsilon_jump=min(1.0d0,exp(met_rat_jump))
END IF
!print*,'epsilon_birth',epsilon_jump,' counter',counter
u=rnunf()
!write(26,100) "birth", iter,nexperts, epsilon_jump,counter
!write(*,100) "birth", iter,nexperts, epsilon_jump,counter
!100 FORMAT(a5,1x,i5,1x,i2,1x,f10.2,1x,i2)
IF(epsilon_jump>u .AND. counter==0) THEN
acc_birth=acc_birth+1
nexperts=nproposed
ndelta=2*(nexperts-1)
DEALLOCATE(beta,delta,tau,prior_ind)
ALLOCATE(beta(nproposed,nbeta),delta(2,nproposed),tau(nproposed),&
prior_ind(nproposed,nchunks))
beta=beta_prop
delta=delta_prop
tau=tau_prop
!CALCULATE NEW PRIOR_IND
DO j=1,nexperts
prior_ind(j,:)=exp(matmul(xdelta,delta(:,j)))/&
sum(exp(matmul(xdelta,delta)),2)
END DO
END IF
DEALLOCATE(tau_prop,beta_prop,delta_prop,prior_ind_prop)
ELSE
!DEATH STEP
nproposed=nexperts-1
ALLOCATE(beta_prop(nproposed,nbeta),delta_prop(2,nproposed),&
prior_ind_prop(nproposed,nchunks),tau_prop(nproposed))
!delta(1,:)=(/0.0, -0.64304866686812/)
!delta(2,:)=(/0.0,0.02764115636602/)
!beta(1,:)=(/3.32283832639074, -73.01996795651158, -19.05586827637581,&
! -34.81691967298258, -46.57188017605456,-8.91738098595595,&
! 59.12663838240293,30.88905387793989,10.01531221750966,&
! -20.12078557915168, -12.60812700514549/)
!beta(2,:)=(/2.91184728926154,66.09989179687619,-26.36129212167623,&
! -13.96764243259529,3.66334777738740,-16.79573458936665,&
! -23.42244397765204,13.10281341818712,22.74223121584767,&
! -8.74595277994963, -8.67254310645845/)
!tau=(/578.022163287414,566.172917759205/)
DO j=1,nproposed
delta_prop(:,j)=delta(:,j)
beta_prop(j,:)=beta(j,:)
tau_prop(j)=tau(j)
END DO
!call wrrrn('beta',beta,nra=nexperts)
!call wrrrn('beta_prop',beta_prop,nra=nproposed)
!call wrrrn('delta',delta,nra=2)
!call wrrrn('delta_prop',delta_prop,nra=2)
!print*,'tau',tau
!print*,'tau_prop',tau_prop
DO j=1,nproposed
prior_ind_prop(j,:)=exp(matmul(xdelta,delta_prop(:,j)))/&
sum(exp(matmul(xdelta,delta_prop)),2)
END DO
DO j=1,nexperts
prior_ind(j,:)=exp(matmul(xdelta,delta(:,j)))/&
sum(exp(matmul(xdelta,delta)),2)
END DO
!Checking dominance condition
counter=0
DO j=1,nproposed
DO k=1,nproposed
IF(maxval(prior_ind_prop(j,:)-prior_ind_prop(k,:))<0.0d0) &
counter=counter+1
END DO
END DO
!print*,'counter',counter
log_like_curr=log_like(nfreq,nchunks,nbeta,nexperts,xx,y,beta,prior_ind)
log_like_new=log_like(nfreq,nchunks,nbeta,nproposed,xx,y,beta_prop,&
prior_ind_prop)
!print*,'log_like_curr', log_like_curr
!print*,'log_like_new', log_like_new
met_rat_jump=log_like_new-log_like_curr
IF(nexperts==2) THEN
epsilon_jump=min(1.0d0,2.0d0*exp(met_rat_jump))
ELSE IF(nexperts==nexp_max) THEN
epsilon_jump=min(1.0d0,0.5d0*exp(met_rat_jump))
ELSE
epsilon_jump=min(1.0d0,exp(met_rat_jump))
END IF
!print*,'epsilon_death',epsilon_jump
u=rnunf()
!write(26,105) "death",iter,nexperts, epsilon_jump, counter
!105 FORMAT(a5,1x,i5,1x,i2,1x,f10.2,1x,i2)
IF(epsilon_jump>u .AND. counter==0) THEN
!IF(epsilon_jump>u) THEN
acc_death=acc_death+1
nexperts=nproposed
ndelta=2*(nexperts-1)
DEALLOCATE(beta,delta,tau,prior_ind)
ALLOCATE(beta(nproposed,nbeta),delta(2,nproposed),tau(nproposed),&
prior_ind(nproposed,nchunks))
beta=beta_prop
delta=delta_prop
tau=tau_prop
!CALCULATE NEW PRIOR_IND
DO j=1,nexperts
prior_ind(j,:)=exp(matmul(xdelta,delta(:,j)))/&
sum(exp(matmul(xdelta,delta)),2)
END DO
END IF
DEALLOCATE(tau_prop,beta_prop,delta_prop,prior_ind_prop)
END IF
!print*,'iteration: ',iter,' nexperts: ',nexperts
IF(nexperts>1) THEN
!Drawing the indicators
ALLOCATE(post_ind(nexperts,nchunks))
call draw_ind(iseed,nfreq,nchunks,nbeta,nexperts,xx,y,beta,&
prior_ind,post_ind,ind)
!do i=1,nexperts
! do j=1,nchunks
! write(15,*) post_ind(i,j)
! end do
!end do
!if(nexperts==4) print*,ind
!IF(mod(iter,100)==0 .AND. iter>nwarmup) THEN
! DO i=1,nchunks
! write(25,*) iter,ind(i)
! END DO
!END IF
DO j=1,nexperts
kk=0
where(ind==j) kk=1
IF(sum(kk)==0) ind(j)=j
END DO
!Drawing the deltas
call draw_delta(iseed,iter,nchunks,ndelta,nexperts,sigdelta,xdelta,&
ind,delta)
END IF
!Drawing the betas
call draw_beta(iseed,nfreq,nchunks,nbeta,nexperts,nbasis,sigalpha,xx,y,&
ind,tau,beta)
call draw_tau(iseed,nbeta,nexperts,nbasis,beta,tau)
!print*,'tau',tau
!prior_ind_hat=0.0
IF(iter>nwarmup) THEN
DO k=1,nchunks
fit_chunk=0.0d0
IF(nexperts>1) THEN
DO j=1,nexperts
fit_chunk=fit_chunk+matmul(xx,beta(j,:))*post_ind(j,k)
!fit_chunk=fit_chunk+matmul(xx,beta(j,:))*prior_ind(j,k)
!prior_ind_hat=prior_ind_hat+prior_ind/dble(nloop-nwarmup)
END DO
fit_ch(:,k)=fit_chunk
ELSE
fit_chunk=matmul(xx,beta(1,:))
fit_ch(:,k)=fit_chunk
END IF
fit(:,k)=fit(:,k)+fit_chunk/dble(nloop-nwarmup)
fit_it(:,iter-nwarmup,k)=fit_chunk
END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!IF(nexperts==3) THEN
!count_fit=count_fit+1
!DO j=1,3
! fit_exp(:,j)=matmul(xx,beta(j,:))
!END DO
!fit_exp_sum=fit_exp_sum+fit_exp
!IF(mod(iter,10)==0) THEN
! DO j=1,nexperts
! DO k=1,nchunks
! write(45,*) post_ind(j,k)
! END DO
! END DO
!END IF
!ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
fit_ch=exp(fit_ch)
sum_fit_ch=sum(fit_ch,1)
DO i=1,nfreq
fit_ch(i,:)=fit_ch(i,:)/sum_fit_ch
END DO
DO k=1,nchunks
cum_fhat(:,k)=cumsum(nfreq,fit_ch(:,k))
END DO
DO i=1,nchunks-1
!dist(i)=sqrt(sum((cum_fhat(:,i)-cum_fhat(:,i+1))**2)/dble(nfreq))
dist(i)=maxval(abs(cum_fhat(:,i)-cum_fhat(:,i+1)))
END DO
!IF(iter>98000) THEN
! DO k=1,nexperts
! DO j=1,nchunks
! write(24,*) prior_ind(k,j)
! END DO
! END DO
!END IF
END IF
IF(nexperts>1) DEALLOCATE(post_ind)
END DO !end main loop
!write(*,*) 'acc_birth: ',acc_birth,' acc_death: ', acc_death
qprop(1)=0.025d0
qprop(2)=0.975d0
DO i=1,nfreq
DO k=1,nchunks
call eqtil(fit_it(i,:,k),nquan,qprop,quan,xlo,xhi,nmiss)
lcl(i,k)=quan(1)
ucl(i,k)=quan(2)
END DO
END DO
DO j=1,nchunks
DO i=1,nfreq
write(20,*) lcl(i,j)
write(22,*) ucl(i,j)
END DO
END DO
!DO k=1,nexperts
! DO j=1,nchunks
! write(24,*) prior_ind_hat(k,j)
! END DO
!END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!fit_exp_sum=fit_exp_sum/dble(count_fit)
!DO j=1,3
! DO i=1,nfreq
! write(40,*) fit_exp_sum(i,j)
! END DO
!END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
DEALLOCATE(tau,beta,delta,prior_ind)
END SUBROUTINE revjump
!
!
FUNCTION log_spect_ar2(nfreq,comp)
IMPLICIT none
INTEGER, INTENT(IN) :: nfreq, comp
DOUBLE PRECISION, DIMENSION(nfreq) :: log_spect_ar2
INTEGER :: i
DOUBLE PRECISION :: phi1,phi2,a,b,c,pi
DOUBLE PRECISION, DIMENSION(nfreq) :: tp
pi=4.0d0*atan(1.0d0)
DO i = 1, nfreq
tp(i)=(dble(i)-1.0d0)/(2.0d0*(dble(nfreq)-1.0d0))
END DO
IF(comp==1) THEN
phi1=0.9d0
phi2=0.0d0
! ELSE IF (comp==2) THEN
! phi1=1.69d0
! phi2=-0.81d0
ELSE
! phi1=1.32d0
! phi2=-0.81d0
phi1=-0.9d0
phi2=0.0d0
END IF
a=1.0d0+phi1**2+phi2**2
b=phi1*phi2-phi1
c=phi2
log_spect_ar2=-log(a+2.0d0*b*cos(2.0d0*pi*tp)-2.0d0*c*cos(4.0d0*pi*tp))
END FUNCTION log_spect_ar2
!
!
FUNCTION cumsum(n,x)
IMPLICIT none
INTEGER, INTENT(IN) :: n
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x
DOUBLE PRECISION, DIMENSION(n) :: cumsum
INTEGER :: i
cumsum(1)=x(1)
DO i=2,n
cumsum(i)=cumsum(i-1)+x(i)
END DO
END FUNCTION cumsum
!
!
SUBROUTINE draw_ind(iseed,nfreq,nchunks,nbeta,nexperts,xx,y,beta,prior_ind,&
post_ind,ind)
IMPLICIT none
INTEGER, INTENT(IN) :: iseed,nfreq,nchunks,nbeta,nexperts
DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: xx
DOUBLE PRECISION, DIMENSION(nfreq,nchunks), INTENT(IN) :: y
DOUBLE PRECISION, DIMENSION(nexperts,nchunks), INTENT(IN) :: prior_ind
DOUBLE PRECISION, DIMENSION(nexperts,nbeta),INTENT(IN) :: beta
INTEGER, DIMENSION(nchunks), INTENT(OUT) :: ind
DOUBLE PRECISION, DIMENSION(nexperts,nchunks),INTENT(OUT) :: post_ind
DOUBLE PRECISION, DIMENSION(nfreq,nexperts) :: ft
DOUBLE PRECISION, DIMENSION(nfreq) :: ydev
DOUBLE PRECISION, DIMENSION(nexperts,nchunks) :: loglike
DOUBLE PRECISION :: u,drnunf
DOUBLE PRECISION :: sum1,sum2
INTEGER :: i,j,k
!Drawing the indicators
DO k=1,nchunks
DO j=1,nexperts
ft(:,j)=matmul(xx,beta(j,:))
ydev=y(:,k)-ft(:,j)
loglike(j,k)=sum(ydev(2:nfreq-1)-exp(ydev(2:nfreq-1)))+&
0.5d0*ydev(1)-0.5d0*exp(ydev(1))+0.5d0*ydev(nfreq)-&
0.5d0*exp(ydev(nfreq))
post_ind(j,k)=loglike(j,k)+log(prior_ind(j,k))
END DO
sum1=sum(exp(post_ind(:,k)))
DO j=1,nexperts
post_ind(j,k)=exp(post_ind(j,k))/sum1
END DO
u=drnunf()
sum2=0.0d0
DO j=1,nexperts
IF(u> sum2 .AND. u<post_ind(j,k)+sum2) THEN
ind(k)=j
END IF
sum2=sum2+post_ind(j,k)
END DO
END DO !end drawing indicators
END SUBROUTINE draw_ind
!
!
SUBROUTINE draw_delta(iseed,iter,nchunks,ndelta,nexperts,sigdelta,xdelta,ind,&
delta)
USE show_int
USE rnunf_int
USE linear_operators
!USE maxdel_m, ONLY:maxdel
IMPLICIT none
INTERFACE
FUNCTION log_target_delta(ndelta,nchunks,nexperts,xdelta,z,delta,sigdelta)
USE show_int
IMPLICIT none
INTEGER, INTENT(IN) :: ndelta,nchunks,nexperts
DOUBLE PRECISION, DIMENSION(nchunks,nexperts), INTENT(IN) :: z
DOUBLE PRECISION, DIMENSION(nchunks,2), INTENT(IN) :: xdelta
DOUBLE PRECISION, DIMENSION(ndelta), INTENT(IN) :: delta
DOUBLE PRECISION, INTENT(IN) :: sigdelta
DOUBLE PRECISION :: log_target_delta
END FUNCTION
FUNCTION log_prop(n,x,mean,var_inv)
IMPLICIT none
INTEGER, INTENT(IN) :: n
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x, mean
DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var_inv
DOUBLE PRECISION :: log_prop
END FUNCTION
FUNCTION mvrnorm(iseed,n,mean,var)
USE lin_eig_self_int
USE linear_operators
IMPLICIT none
INTEGER, INTENT(IN) :: n,iseed
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: mean
DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var
DOUBLE PRECISION, DIMENSION(n) :: mvrnorm
END FUNCTION
FUNCTION mvrnorm1(iseed,n,mean,var)
USE linear_operators
USE show_int
IMPLICIT none
INTEGER, INTENT(IN) :: n,iseed
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: mean
DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var
DOUBLE PRECISION, DIMENSION(n) :: mvrnorm1
END FUNCTION
FUNCTION mvrnorm2(iseed,n,mean,var) RESULT(r)
!n is the number of draws from the multivariate normal
USE rnmvn_int
USE chfac_int
IMPLICIT none
INTEGER, PARAMETER :: DP = kind(1.d0)
INTEGER :: iseed, n
REAL(DP) :: mean(:), var(:,:)
REAL (DP) :: rsig(size(mean),size(mean)), r(n,size(mean))
INTEGER :: irank
END FUNCTION
END INTERFACE
INTEGER, INTENT(IN) :: nchunks,ndelta,nexperts,iseed
DOUBLE PRECISION, INTENT(IN) :: sigdelta
DOUBLE PRECISION, DIMENSION(nchunks,2), INTENT(IN) :: xdelta
INTEGER, DIMENSION(nchunks), INTENT(IN) :: ind
INTEGER, INTENT(IN) :: iter
DOUBLE PRECISION, DIMENSION(2,nexperts), INTENT(INOUT) :: delta
DOUBLE PRECISION, DIMENSION(2,ndelta) :: temp_delta
DOUBLE PRECISION, DIMENSION(ndelta) :: temp1
DOUBLE PRECISION, DIMENSION(ndelta) :: delta_mean
DOUBLE PRECISION, DIMENSION(ndelta,ndelta) :: delta_var
DOUBLE PRECISION, DIMENSION(ndelta) :: sold,sold1
DOUBLE PRECISION, DIMENSION(ndelta) :: derd
DOUBLE PRECISION, DIMENSION(ndelta,ndelta) :: hessd
DOUBLE PRECISION, DIMENSION(nchunks,nexperts) :: z
DOUBLE PRECISION :: met_rat_delta
INTEGER :: reject
DOUBLE PRECISION :: tolf,tolx,tol,dmach,u
DOUBLE PRECISION :: log_target_delta_old,log_target_delta_new
!DOUBLE PRECISION :: drnunf,epsilond
DOUBLE PRECISION :: epsilond
DOUBLE PRECISION :: log_prop_delta_old,log_prop_delta_new
DOUBLE PRECISION, DIMENSION(ndelta) :: pard
INTEGER, DIMENSION(nchunks) :: kk
INTEGER :: irank,nr
INTEGER :: i,j,k,jj
INTEGER :: count, counter
DOUBLE PRECISION :: rnorm(1,ndelta),prior_ind_prop(nexperts,nchunks)
DOUBLE PRECISION :: delta_prop(2,nexperts)
!Drawing the deltas
z=0.0d0
DO i=1,nchunks
DO j=1,nexperts
IF(ind(i)==j) z(i,j)=1.0d0
END DO
END DO
count=0
DO j=2,nexperts
DO k=1,2
count=count+1
temp_delta(1,count)=delta(k,j)
END DO
END DO
!IF(iter<100) THEN
!call drnun(ndelta,pard)
!call drnnor(ndelta,pard)
!pard=0.1d0*pard
pard=0.0d0
!ELSE
! pard=temp_delta(1,:)
!END IF
!call maxdel(ndelta,pard,ind,nchunks,nexperts,xdelta,sold1)
!print*,'opt',sold1
tolf=1.0d0
tolx=1.0d0
jj=0
DO WHILE(tolf>.01d0 .AND. tolx>.01d0 .AND. jj<=20)
jj=jj+1
call derivsd(z,ndelta,nexperts,nchunks,pard,xdelta,sigdelta,derd,hessd)
call dlslsf(ndelta,hessd,ndelta,derd,sold)
pard=pard-sold
tolf=sum(abs(derd))
tolx=sum(abs(-hessd.ix.derd))
END DO
delta_mean=pard
!print*,'nr',delta_mean
delta_var=-.i.hessd
rnorm=mvrnorm2(iseed,1,delta_mean,delta_var)
temp_delta(2,:)=rnorm(1,:)
!temp_delta(2,:)=mvrnorm1(iseed,ndelta,delta_mean,delta_var)
log_prop_delta_old=log_prop(ndelta,temp_delta(1,:),delta_mean,-hessd)
log_prop_delta_new=log_prop(ndelta,temp_delta(2,:),delta_mean,-hessd)
log_target_delta_old=log_target_delta(ndelta,nchunks,nexperts,xdelta,z,&
temp_delta(1,:),sigdelta)
!print*,'log_target_delta_old',log_target_delta_old
log_target_delta_new=log_target_delta(ndelta,nchunks,nexperts,xdelta,z,&
temp_delta(2,:),sigdelta)
met_rat_delta=log_target_delta_new-log_target_delta_old+log_prop_delta_old-&
log_prop_delta_new
!
count=0
delta_prop=0.0
DO j=2,nexperts
DO k=1,2
count=count+1
delta_prop(k,j)=temp_delta(2,count)
END DO
END DO
DO j=1,nexperts
prior_ind_prop(j,:)=exp(matmul(xdelta,delta_prop(:,j)))/&
sum(exp(matmul(xdelta,delta_prop)),2)
END DO
!Checking dominance condition
counter=0
DO j=1,nexperts
DO k=1,nexperts
IF(maxval(prior_ind_prop(j,:)-prior_ind_prop(k,:))<0.0d0) &
counter=counter+1
END DO
END DO
u=rnunf()
epsilond=min(1.0d0,exp(met_rat_delta))
!print*,'epsilond',epsilond,'counter',counter
!write(26,105) "gating",epsilond,counter
!105 FORMAT(a5,1x,f10.2,1x,i2)
IF(epsilond > u .AND. counter==0) THEN
count=0
DO j=2,nexperts
DO k=1,2
count=count+1
delta(k,j)=temp_delta(2,count)
END DO
END DO
END IF
!IF(u>epsilond) THEN
! temp_delta(2,:)=temp_delta(1,:)
!END IF
END SUBROUTINE draw_delta
!
!
SUBROUTINE draw_beta(iseed,nfreq,nchunks,nbeta,nexperts,nbasis,sigalpha,xx,y,&
ind,tau,beta)
USE show_int
USE wrrrn_int
USE rnunf_int
USE lin_eig_self_int
USE linear_operators
USE maxb_m, ONLY: maxb
IMPLICIT none
INTERFACE
FUNCTION findloc(m,n,ind,ii)
INTEGER, INTENT(IN) :: m,n
INTEGER, DIMENSION(n), INTENT(IN) :: ind
INTEGER, DIMENSION(m) :: findloc
INTEGER :: ii,j,k
END FUNCTION
FUNCTION log_prop(n,x,mean,var_inv)
IMPLICIT none
INTEGER, INTENT(IN) :: n
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x, mean
DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var_inv
DOUBLE PRECISION :: log_prop
END FUNCTION
FUNCTION log_target(nfreq,ncoly,nbeta,x,y,par,sigalpha,tau)
IMPLICIT none
INTEGER, INTENT(IN) :: nfreq,ncoly,nbeta
DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: x
DOUBLE PRECISION, DIMENSION(nfreq,ncoly), INTENT(IN) :: y
DOUBLE PRECISION, DIMENSION(nbeta), INTENT(IN) :: par
DOUBLE PRECISION, INTENT(IN) :: sigalpha, tau
DOUBLE PRECISION :: log_target
END FUNCTION
FUNCTION mvrnorm1(iseed,n,mean,var)
USE linear_operators
IMPLICIT none
INTEGER, INTENT(IN) :: n,iseed
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: mean
DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var
DOUBLE PRECISION, DIMENSION(n) :: mvrnorm1
END FUNCTION
FUNCTION mvrnorm(iseed,n,mean,var)
USE lin_eig_self_int
USE linear_operators
IMPLICIT none
INTEGER, INTENT(IN) :: n,iseed
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: mean
DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var
DOUBLE PRECISION, DIMENSION(n) :: mvrnorm
END FUNCTION
FUNCTION mvrnorm2(iseed,n,mean,var) RESULT(r)
!n is the number of draws from the multivariate normal
USE rnmvn_int
USE chfac_int
IMPLICIT none
INTEGER, PARAMETER :: DP = kind(1.d0)
INTEGER :: iseed, n
REAL(DP) :: mean(:), var(:,:)
REAL (DP) :: rsig(size(mean),size(mean)), r(n,size(mean))
INTEGER :: irank
END FUNCTION
END INTERFACE
INTEGER, INTENT(IN) :: nfreq,nchunks,nbeta,nexperts,nbasis,iseed
DOUBLE PRECISION, INTENT(IN) :: sigalpha
DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: xx
DOUBLE PRECISION, DIMENSION(nfreq,nchunks), INTENT(IN) :: y
DOUBLE PRECISION, DIMENSION(nexperts),INTENT(IN) :: tau
INTEGER, DIMENSION(nchunks),INTENT(IN) :: ind
DOUBLE PRECISION, DIMENSION(nexperts,nbeta), INTENT(INOUT) :: beta
DOUBLE PRECISION, DIMENSION(nexperts,nbeta,2) :: temp_beta
DOUBLE PRECISION, DIMENSION(nfreq,nchunks) :: fit
DOUBLE PRECISION, DIMENSION(nfreq) :: fit_chunk
DOUBLE PRECISION, DIMENSION(nbeta) :: derb
DOUBLE PRECISION, DIMENSION(nbeta) :: param
DOUBLE PRECISION, DIMENSION(nbeta) :: temp
DOUBLE PRECISION, DIMENSION(nbeta,nbeta) :: hessb_inv
DOUBLE PRECISION, DIMENSION(nbeta) :: beta_mean
DOUBLE PRECISION, DIMENSION(nbeta,nbeta) :: beta_var
DOUBLE PRECISION, DIMENSION(nbeta) :: solb
DOUBLE PRECISION, DIMENSION(nbeta) :: sol
DOUBLE PRECISION, DIMENSION(nbeta) :: rnor
DOUBLE PRECISION, DIMENSION(nbeta,nbeta) :: hessb,beta_var_inv
DOUBLE PRECISION, DIMENSION(nbeta,nbeta) :: rsig
DOUBLE PRECISION, DIMENSION(nfreq,nexperts) :: ft_hat,ft
DOUBLE PRECISION, DIMENSION(nfreq,nchunks) :: fhat
DOUBLE PRECISION, DIMENSION(nfreq) :: ydev
DOUBLE PRECISION, DIMENSION(nexperts,nchunks) :: loglike,post_ind,&
prior_ind
DOUBLE PRECISION, DIMENSION(nexperts) :: met_rat
INTEGER :: reject
DOUBLE PRECISION, DIMENSION(nexperts,nbeta) :: beta0
DOUBLE PRECISION :: tolf,tolx,tol,dmach,log_prop_old,log_prop_new,u
!DOUBLE PRECISION :: log_target_old,log_target_new,epsilonb,drnunf
DOUBLE PRECISION :: log_target_old,log_target_new,epsilonb
INTEGER, DIMENSION(nchunks) :: kk
INTEGER, DIMENSION(nbeta) :: nev
INTEGER, ALLOCATABLE :: kk_ind(:)
DOUBLE PRECISION, ALLOCATABLE :: hh(:)
INTEGER :: i,j,k,iter,ii,jj
INTEGER :: n_chunk_ind,n_exp
DOUBLE PRECISION, DIMENSION(nbeta,nbeta) :: eigenvec
DOUBLE PRECISION, DIMENSION(nbeta) :: eigenval
DOUBLE PRECISION :: rnorm(1,nbeta)
temp_beta(:,:,1)=beta
DO j=1,nexperts
kk=0
where(ind.eq.j) kk=1
n_chunk_ind=sum(kk)
ALLOCATE(kk_ind(n_chunk_ind))
!find the chunks belonging to expert j
kk_ind=findloc(n_chunk_ind,nchunks,ind,j)
!generate beta
!find the posterior mode of beta
!param=temp_beta(j,:,1)
param=0.0d0
tolf=1.0d0
tolx=1.0d0
ii=0
DO WHILE(tolf>.001d0 .AND. tolx>.00001d0)
ii=ii+1
nev=0
call derivs(nbeta,nfreq,n_chunk_ind,param,xx,y(:,kk_ind),&
sigalpha,tau(j),derb,hessb)
call lin_eig_self(hessb,eigenval,v=eigenvec)
WHERE(eigenval.LE.0.0d0) nev=1
IF(SUM(nev)>0) EXIT
call dlslsf(nbeta,hessb,nbeta,derb,solb)
param=param-solb
tolf=sum(abs(derb))
tolx=sum(abs(-hessb.ix.derb))
END DO
IF(SUM(nev)>0) THEN
param=temp_beta(j,:,1)
ALLOCATE(hh(n_chunk_ind))
hh=1.0d0
call maxb(param,tau(j),hh,sigalpha,nfreq,n_chunk_ind,nbeta,&
nexperts,nbasis,xx,y(:,kk_ind),sol)
call derivs(nbeta,nfreq,n_chunk_ind,sol,xx,y(:,kk_ind),&
sigalpha,tau(j),derb,hessb)
DEALLOCATE(hh)
END IF
beta_mean=param
beta_var_inv=-hessb/1.2d0
beta_var=-1.2d0*.i.hessb
!temp_beta(j,:,2)=mvrnorm1(iseed,nbeta,beta_mean,beta_var)
rnorm=mvrnorm2(iseed,1,beta_mean,beta_var)
!print*,'rnorm',rnorm
temp_beta(j,:,2)=rnorm(1,:)
log_prop_old=log_prop(nbeta,temp_beta(j,:,1),beta_mean,beta_var_inv)
!beta(1,:)=(/2.95326589079431, -52.28456352975171, -59.68971488606051,&
!-14.03614907332720, -15.24975538988489, -39.33304649805422,&
! 47.47076878242155, 43.85075118448089, 16.11726555949620,&
! -15.46952991800888, -11.25003783572547/)
log_prop_new=log_prop(nbeta,temp_beta(j,:,2),beta_mean,beta_var_inv)
log_target_old=log_target(nfreq,n_chunk_ind,nbeta,xx,y(:,kk_ind),&
temp_beta(j,:,1),sigalpha,tau(j))
log_target_new=log_target(nfreq,n_chunk_ind,nbeta,xx,y(:,kk_ind),&
temp_beta(j,:,2),sigalpha,tau(j))
!Metropolis ratio
met_rat(j)=log_target_new-log_target_old+log_prop_old-log_prop_new
u=rnunf()
epsilonb=min(1.0d0,exp(met_rat(j)))
IF(u>epsilonb) THEN
temp_beta(j,:,2)=temp_beta(j,:,1)
!IF(j==2) reject=reject+1
END IF
beta(j,:)=temp_beta(j,:,2)
DEALLOCATE(kk_ind)
END DO !end expert loop
END SUBROUTINE draw_beta
!
!
SUBROUTINE draw_tau(iseed,nbeta,nexperts,nbasis,beta,tau)
USE show_int
USE gamin_int
USE gamdf_int
USE rngam_int
USE rnunf_int
IMPLICIT none
INTEGER, INTENT(IN) :: nbeta,nexperts,nbasis,iseed
DOUBLE PRECISION, DIMENSION(nexperts,nbeta), INTENT(IN) :: beta
DOUBLE PRECISION, DIMENSION(nexperts),INTENT(INOUT) :: tau
!DOUBLE PRECISION :: taua,taub,gamdev,const1,pinv,dgamin,dgamdf
DOUBLE PRECISION :: taua,taub,gamdev(1),const1,pinv
!DOUBLE PRECISION :: u,drnunf
DOUBLE PRECISION :: u
INTEGER :: j
!Drawing tau
DO j=1,nexperts
taua=0.5d0*dble(nbasis)
taub=0.5d0*sum(beta(j,2:nbeta)**2)
IF (j==1) THEN
!call drngam(1,taua,gamdev)
call rngam(taua,gamdev)
tau(j)=taub/gamdev(1)
!Drawing Tauj constrained to be less than tauj-1
ELSE
!const1=dgamdf(taub/tau(j-1),taua)
const1=gamdf(taub/tau(j-1),taua)
!u=drnunf()
u=rnunf()
!pinv=u*(1.0d0-const1)+const1
!tau(j)=taub/dgamin(pinv,taua)
pinv=1.0d0-u*(1.0d0-const1)
IF(pinv>0.999d0) THEN
tau(j)=tau(j-1)-1.0d0
ELSE
!tau(j)=taub/dgamin(pinv,taua)
tau(j)=taub/gamin(pinv,taua)
END IF
END IF
END DO
END SUBROUTINE draw_tau
!
!
SUBROUTINE em(iseed,nfreq,nchunks,nbeta,ndelta,nexperts,nbasis,sigalpha,&
sigdelta,xx,xdelta,y,ind0,beta0,prior_ind)
USE show_int
USE maxb_m, ONLY: maxb
USE maxd_m, ONLY: maxd
IMPLICIT none
INTEGER, INTENT(IN) :: iseed,nfreq,nchunks,nbeta,ndelta,nexperts,nbasis
DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: xx
DOUBLE PRECISION, DIMENSION(nchunks,2), INTENT(IN) :: xdelta
DOUBLE PRECISION, DIMENSION(nfreq,nchunks), INTENT(IN) :: y
DOUBLE PRECISION, INTENT(IN) :: sigalpha,sigdelta
DOUBLE PRECISION, DIMENSION(nexperts,nchunks) :: h
DOUBLE PRECISION, DIMENSION(nchunks) :: h1
INTEGER, DIMENSION(nchunks), INTENT(OUT) :: ind0
DOUBLE PRECISION, DIMENSION(nexperts,nbeta), INTENT(OUT) :: beta0
DOUBLE PRECISION, DIMENSION(nexperts,nchunks), INTENT(OUT) :: prior_ind
DOUBLE PRECISION, ALLOCATABLE :: tau(:),beta_curr(:,:),delta(:,:)
DOUBLE PRECISION, ALLOCATABLE :: loglike(:,:)
DOUBLE PRECISION, ALLOCATABLE :: ydev(:),ft(:,:),denom(:),&
h_sum(:),param(:), sol(:),par(:),sold(:)
INTEGER :: i,j,k,mm,iter,jj,kk
DOUBLE PRECISION :: f,tolf,fd
DOUBLE PRECISION :: crit
ALLOCATE(beta_curr(2,nbeta*nexperts),ft(nfreq,nexperts),ydev(nfreq), &
loglike(nexperts,nchunks),denom(nchunks),delta(2,nexperts),&
h_sum(nchunks),tau(nexperts),param(nbeta),sol(nbeta),par(ndelta),&
sold(ndelta))
IF(nexperts>=2) THEN
beta_curr=0.0d0
tau(1)=10000.0d0
DO j=2,nexperts
tau(j)=tau(j-1)/10.0d0
END DO
delta=0.0d0
iter=0
crit=100.0d0
!DO mm=1,5 !start EM loop
DO WHILE(crit>.00001d0 .AND. iter <=20) !start EM loop
iter=iter+1
!E-step
DO k=1,nchunks
DO j=1,nexperts
ft(:,j)=matmul(xx,beta_curr(1,1+(j-1)*nbeta:j*nbeta))
ydev=y(:,k)-ft(:,j)
loglike(j,k)=sum(ydev(2:nfreq-1)-exp(ydev(2:nfreq-1)))+&
0.5*(ydev(1)-exp(ydev(1)))+&
0.5*(ydev(nfreq)-exp(ydev(nfreq)))
END DO
END DO
denom=0.0d0
DO j=1,nexperts
denom=denom+exp(matmul(xdelta,delta(:,j)))
END DO
DO j=1,nexperts
prior_ind(j,:)=exp(matmul(xdelta,delta(:,j)))/denom
END DO
h=prior_ind*exp(loglike)
h_sum=sum(h,1)
DO j=1,nexperts
h(j,:)=h(j,:)/h_sum
END DO
IF(iter==1) THEN
DO j=1,nexperts
call drnun(nchunks,h(j,:))
END DO
h_sum=sum(h,1)
DO j=1,nexperts
h(j,:)=h(j,:)/h_sum
END DO
END IF
!M-step
!wrt the betas
DO j=1,nexperts
param=beta_curr(1,1+(j-1)*nbeta:j*nbeta)
call maxb(param,tau(j),h(j,:),sigalpha,nfreq,nchunks,nbeta,&
nexperts,nbasis,xx,y,sol)
beta_curr(2,1+(j-1)*nbeta:j*nbeta)=sol
END DO
!call show(beta_curr)
!wrt the deltas
par=0.0d0
call maxd(ndelta,par,h,sigdelta,nchunks,nexperts,xdelta,sold)
! print*,sold
DO j=1,nexperts-1
delta(:,j+1)=sold(1+(j-1)*2:j*2)
END DO
!wrt the taus
DO j=1,nexperts
tau(j)=0.5d0*sum(beta_curr(2,2+(j-1)*nbeta:j*nbeta)**2)
END DO
!print*,tau
crit=maxval(abs(beta_curr(1,:)-beta_curr(2,:))/(beta_curr(1,:)+.001d0))
beta_curr(1,:)=beta_curr(2,:)
END DO !end of EM loop
DO j=1,nexperts
beta0(j,:)=beta_curr(1,1+(j-1)*nbeta:j*nbeta)
END DO
ind0=maxloc(h,dim=1)
ELSE
h1=1.0d0
call maxb(param,tau(j),h1,sigalpha,nfreq,nchunks,nbeta,&
nexperts,nbasis,xx,y,sol)
beta0(1,:)=sol
ind0=1
prior_ind(1,:)=1.0d0
END IF
DEALLOCATE(beta_curr,ft,loglike,ydev,denom,delta,h_sum,param,&
tau,sol,par,sold)
END SUBROUTINE em
!
!
SUBROUTINE funcb(nbeta,param,f)
USE maxb_m, ONLY:tau_exp,pp,sigalpha,nfreq,nchunks,nexperts,nbasis,xx,y
IMPLICIT none
INTEGER, INTENT(IN) :: nbeta
DOUBLE PRECISION, DIMENSION(nbeta), INTENT(IN) :: param
DOUBLE PRECISION, INTENT(OUT) :: f
DOUBLE PRECISION, ALLOCATABLE :: ydev(:,:),xxparam(:)
INTEGER :: i,j,k
ALLOCATE(xxparam(nfreq),ydev(nfreq,nchunks))
xxparam=matmul(xx,param)
DO j=1,nchunks
ydev(:,j)=y(:,j)-xxparam
END DO
f=0.0d0
DO k=1,nchunks
f=f+pp(k)*(0.5*ydev(1,k)-0.5*exp(ydev(1,k)) &
+0.5*ydev(nfreq,k)-0.5*exp(ydev(nfreq,k)) &
+sum(ydev(2:nfreq-1,k)-exp(ydev(2:nfreq-1,k))))
END DO
f=f-0.5*(sum(param(2:nbeta)**2)/tau_exp+param(1)**2/sigalpha);
f=-f
DEALLOCATE(xxparam,ydev)
END SUBROUTINE funcb
!
!
SUBROUTINE gradb(nbeta,param,g)
USE maxb_m, ONLY:tau_exp,pp,sigalpha,nfreq,nchunks,nexperts,xx,y
IMPLICIT none
INTEGER, INTENT(IN) :: nbeta
DOUBLE PRECISION, DIMENSION(nbeta), INTENT(IN) :: param
DOUBLE PRECISION, DIMENSION(nbeta),INTENT(OUT) :: g
INTEGER :: i,j,k
g=0.0d0
g(1)=-param(1)/sigalpha
g(2:nbeta)=-param(2:nbeta)/tau_exp
DO i=1,nfreq
DO k=1,nchunks
IF(i>1 .AND. i<nfreq) THEN
g=g-pp(k)*xx(i,:)+pp(k)*xx(i,:)*exp(y(i,k)-sum(xx(i,:)*param))
ELSE
g=g-0.5d0*pp(k)*xx(i,:)+0.5*pp(k)*xx(i,:) &
*exp(y(i,k)-sum(xx(i,:)*param))
END IF
END DO
END DO
g=-g
END SUBROUTINE gradb
!
!
SUBROUTINE hessb(nbeta,param,hess,ldh)
USE maxb_m, ONLY:tau_exp,pp,sigalpha,nfreq,nchunks,nexperts,nbasis,xx,y
IMPLICIT none
INTERFACE
FUNCTION out_prod(n,x)
INTEGER, INTENT(IN) :: n
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x
DOUBLE PRECISION, DIMENSION(n,n) :: out_prod
END FUNCTION
END INTERFACE
INTEGER, INTENT(IN) :: nbeta,ldh
DOUBLE PRECISION, DIMENSION(nbeta), INTENT(IN) :: param
DOUBLE PRECISION, DIMENSION(nbeta,nbeta),INTENT(OUT) :: hess
INTEGER :: i,j,k
hess=0.0d0
hess(1,1)=-1.0/sigalpha
hess(2:nbeta,2:nbeta)=-1.0d0/tau_exp
DO i=1,nfreq
DO k=1,nchunks
IF(i>1 .AND. i<nfreq) THEN
hess=hess-pp(k)*out_prod(nbeta,xx(i,:))*&
exp(y(i,k)-sum(xx(i,:)*param))
ELSE
hess=hess-pp(k)*out_prod(nbeta,xx(i,:))*&
0.5*exp(y(i,k)-sum(xx(i,:)*param))
END IF
END DO
END DO
hess=-hess
END SUBROUTINE hessb
!
!
FUNCTION out_prod(n,x)
IMPLICIT none
INTEGER, INTENT(IN) :: n
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x
DOUBLE PRECISION, DIMENSION(n,n) :: out_prod
INTEGER :: i
DO i=1,n
out_prod(i,:)=x(i)*x
END DO
END FUNCTION out_prod
!
!
SUBROUTINE funcd(ndelta,param,f)
USE maxd_m, ONLY:h,sigdelta,nchunks,nexperts,xdelta
IMPLICIT none
INTEGER, INTENT(IN) :: ndelta
DOUBLE PRECISION, DIMENSION(ndelta), INTENT(IN) :: param
DOUBLE PRECISION, INTENT(OUT) :: f
DOUBLE PRECISION, DIMENSION(nexperts*2) :: par1
DOUBLE PRECISION, DIMENSION(nchunks) :: denom
DOUBLE PRECISION, DIMENSION(nexperts,nchunks) :: prior_ind
INTEGER :: i,j,k
par1=0.0d0
par1(3:2*nexperts)=param
denom=0.0d0
DO j=1,nexperts
denom=denom+exp(matmul(xdelta,par1(1+(j-1)*2:2*j)))
END DO
DO j=1,nexperts
prior_ind(j,:)=exp(matmul(xdelta,par1(1+(j-1)*2:2*j)))/denom
END DO
f=SUM(h*log(prior_ind))-0.5d0*SUM(param**2)/sigdelta
f=-f
END SUBROUTINE funcd
!
!
SUBROUTINE gradd(ndelta,param,g)
USE maxd_m, ONLY:h,sigdelta,nchunks,nexperts,xdelta
IMPLICIT none
INTEGER, INTENT(IN) :: ndelta
DOUBLE PRECISION, DIMENSION(ndelta), INTENT(IN) :: param
DOUBLE PRECISION, DIMENSION(ndelta),INTENT(OUT) :: g
DOUBLE PRECISION, DIMENSION(nexperts*2) :: par1
DOUBLE PRECISION, DIMENSION(nchunks) :: denom
DOUBLE PRECISION, DIMENSION(nexperts,nchunks) :: prior_ind
INTEGER :: i,j,k
par1=0.0d0
par1(3:2*nexperts)=param
denom=0.0d0
DO j=1,nexperts
denom=denom+exp(matmul(xdelta,par1(1+(j-1)*2:2*j)))
END DO
DO j=1,nexperts
prior_ind(j,:)=exp(matmul(xdelta,par1(1+(j-1)*2:2*j)))/denom
END DO
g=-par1(3:2*nexperts)/sigdelta
DO k=1,nexperts-1
DO i=1,nchunks
g(1+(k-1)*2:k*2)=g(1+(k-1)*2:k*2)+(h(k,i)-prior_ind(k,i))*xdelta(i,:)
END DO
END DO
g=-g
END SUBROUTINE gradd
!
!
SUBROUTINE hessd(ndelta,param,hd,ldh)
USE maxd_m, ONLY:h,sigdelta,nchunks,nexperts,xdelta
IMPLICIT none
INTERFACE
FUNCTION out_prod(n,x)
INTEGER, INTENT(IN) :: n
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x
DOUBLE PRECISION, DIMENSION(n,n) :: out_prod
END FUNCTION
END INTERFACE
INTEGER, INTENT(IN) :: ndelta,ldh
DOUBLE PRECISION, DIMENSION(ndelta), INTENT(IN) :: param
DOUBLE PRECISION, DIMENSION(ndelta,ndelta),INTENT(OUT) :: hd
DOUBLE PRECISION, DIMENSION(nexperts*2) :: par1
DOUBLE PRECISION, DIMENSION(nchunks) :: denom
DOUBLE PRECISION, DIMENSION(nexperts,nchunks) :: prior_ind
INTEGER :: i,j,k,l,kron
par1=0.0d0
par1(3:2*nexperts)=param
denom=0.0d0
DO j=1,nexperts
denom=denom+exp(matmul(xdelta,par1(1+(j-1)*2:2*j)))
END DO
DO j=1,nexperts
prior_ind(j,:)=exp(matmul(xdelta,par1(1+(j-1)*2:2*j)))/denom
END DO
hd=0.0d0
DO i=1,(nexperts-1)*2
hd(i,i)=-1.0d0/sigdelta
END DO
DO k=1,nexperts-1
DO l=1,nexperts-1
IF(k==l) THEN
kron=1
ELSE
kron=0
END IF
DO i=1,nchunks
hd(1+(k-1)*2:2*k,1+(l-1)*2:2*l)=hd(1+(k-1)*2:2*k,1+(l-1)*2:2*l)-&
prior_ind(k,i)*(kron-prior_ind(l,i))*out_prod(2,xdelta(i,:))
END DO
END DO
END DO
hd=-hd
END SUBROUTINE hessd
!
!
FUNCTION findloc(m,n,ind,ii)
!find the locations of all the values of
!ind equal to ii
IMPLICIT none
INTEGER, INTENT(IN) :: m,n
INTEGER, DIMENSION(n), INTENT(IN) :: ind
INTEGER, DIMENSION(m) :: findloc
INTEGER :: ii,j,k
k=1
DO j=1,n
IF(ind(j).eq. ii) THEN
findloc(k)=j
k=k+1
END IF
END DO
END FUNCTION findloc
!
!
SUBROUTINE derivs(nbeta,nfreq,ncoly,beta,x,y,sigalpha,tau,der,hess)
USE show_int
IMPLICIT none
INTERFACE
FUNCTION out_prod(n,x)
INTEGER, INTENT(IN) :: n
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x
DOUBLE PRECISION, DIMENSION(n,n) :: out_prod
END FUNCTION
END INTERFACE
INTEGER, INTENT(IN) :: nbeta,nfreq,ncoly
DOUBLE PRECISION, DIMENSION(nbeta), INTENT(IN) :: beta
DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: x
DOUBLE PRECISION, DIMENSION(nfreq,ncoly), INTENT(IN) :: y
DOUBLE PRECISION, INTENT(IN) :: sigalpha,tau
DOUBLE PRECISION, DIMENSION(nbeta), INTENT(OUT) :: der
DOUBLE PRECISION, DIMENSION(nbeta,nbeta), INTENT(OUT) :: hess
INTEGER :: i,k
DOUBLE PRECISION :: exp_ys
der=0.0d0
der(1)=-beta(1)/sigalpha
der(2:nbeta)=-beta(2:nbeta)/tau
hess=0.0d0
hess(1,1)=-1.0d0/sigalpha
DO i=2,nbeta
hess(i,i)=-1.0d0/tau
END DO
DO i=1,nfreq
DO k=1,ncoly
!exp_ys=exp(y(i,k)-sum(x(i,:)*beta))
IF(i>1 .and. i<nfreq) THEN
der=der-x(i,:)+x(i,:)*exp(y(i,k)-sum(x(i,:)*beta))
hess=hess-out_prod(nbeta,x(i,:))*exp(y(i,k)-sum(x(i,:)*beta))
ELSE
der=der-0.5d0*x(i,:)+0.5d0*x(i,:)*exp(y(i,k)-sum(x(i,:)*beta))
hess=hess-0.5d0*out_prod(nbeta,x(i,:))*&
exp(y(i,k)-sum(x(i,:)*beta))
END IF
END DO
END DO
END SUBROUTINE derivs
!
!
SUBROUTINE derivs1(nbeta,nfreq,ncoly,beta,x,y,sigalpha,tau,der,hess)
USE show_int
IMPLICIT none
INTERFACE
FUNCTION out_prod(n,x)
INTEGER, INTENT(IN) :: n
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x
DOUBLE PRECISION, DIMENSION(n,n) :: out_prod
END FUNCTION
END INTERFACE
INTEGER, INTENT(IN) :: nbeta,nfreq,ncoly
DOUBLE PRECISION, DIMENSION(nbeta), INTENT(IN) :: beta
DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: x
DOUBLE PRECISION, DIMENSION(nfreq,ncoly), INTENT(IN) :: y
DOUBLE PRECISION, INTENT(IN) :: sigalpha,tau
DOUBLE PRECISION, DIMENSION(nbeta), INTENT(OUT) :: der
DOUBLE PRECISION, DIMENSION(nbeta,nbeta), INTENT(OUT) :: hess
INTEGER :: i,k
DOUBLE PRECISION :: exp_ys
der=0.0d0
der(1)=-beta(1)/sigalpha
der(2:nbeta)=-beta(2:nbeta)/tau
hess=0.0d0
hess(1,1)=-1.0d0/sigalpha
DO i=2,nbeta
hess(i,i)=-1.0d0/tau
END DO
call show(y)
DO i=1,nfreq
DO k=1,ncoly
IF (mod(i,nfreq)==0.OR.mod(i,nfreq)==1) THEN
der=der-0.5d0*x(i,:)+0.5d0*x(i,:)*exp(y(i,k)-sum(x(i,:)*beta))
hess=hess-out_prod(nbeta,x(i,:))*0.5d0*&
exp(y(i,k)-sum(x(i,:)*beta))
ELSE
der=der-x(i,:)+x(i,:)*exp(y(i,k)-sum(x(i,:)*beta))
hess=hess-out_prod(nbeta,x(i,:))*exp(y(i,k)-sum(x(i,:)*beta))
END IF
END DO
END DO
END SUBROUTINE derivs1
!
!
SUBROUTINE derivsd(z,ndelta,nexperts,nchunks,param,xdelta,sigdelta,der,hess)
USE show_int
USE linear_operators
IMPLICIT none
INTERFACE
FUNCTION out_prod(n,x)
INTEGER, INTENT(IN) :: n
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x
DOUBLE PRECISION, DIMENSION(n,n) :: out_prod
END FUNCTION
END INTERFACE
INTEGER, INTENT(IN) :: nexperts,nchunks,ndelta
!INTEGER, INTENT(IN),DIMENSION(nchunks) :: ind
DOUBLE PRECISION, INTENT(IN),DIMENSION(nchunks,nexperts) :: z
DOUBLE PRECISION, DIMENSION(nchunks,2), INTENT(IN) :: xdelta
DOUBLE PRECISION, DIMENSION(ndelta), INTENT(IN) :: param
DOUBLE PRECISION, INTENT(IN) :: sigdelta
DOUBLE PRECISION, DIMENSION(ndelta), INTENT(OUT) :: der
DOUBLE PRECISION, DIMENSION(ndelta,ndelta), INTENT(OUT) :: hess
DOUBLE PRECISION, DIMENSION(nexperts,nchunks) :: prior_ind
DOUBLE PRECISION, DIMENSION(nchunks) :: denom
DOUBLE PRECISION, DIMENSION(nexperts*2) :: param1
DOUBLE PRECISION :: kron
INTEGER :: i,j,k,l
param1=0.0d0
param1(3:nexperts*2)=param
denom=0.0d0
DO j=1,nexperts
denom=denom+exp(matmul(xdelta,param1(1+(j-1)*2:2*j)))
END DO
DO j=1,nexperts
prior_ind(j,:)=exp(matmul(xdelta,param1(1+(j-1)*2:2*j)))/denom
END DO
der=-param/sigdelta
hess=-1.0d0/sigdelta*eye(ndelta)
DO k=2,nexperts
DO i=1,nchunks
der(1+(k-2)*2:2*(k-1))=der(1+(k-2)*2:2*(k-1))+(z(i,k)-prior_ind(k,i))*&
xdelta(i,:)
END DO
END DO
DO k=2,nexperts
DO l=2,nexperts
IF(k==l) THEN
kron=1.0d0
ELSE
kron=0.0d0
END IF
DO i=1,nchunks
hess(1+(k-2)*2:2*(k-1),1+(l-2)*2:2*(l-1))=&
hess(1+(k-2)*2:2*(k-1),1+(l-2)*2:2*(l-1))-prior_ind(k,i)*&
(kron-prior_ind(l,i))*out_prod(2,xdelta(i,:))
END DO
END DO
END DO
END SUBROUTINE derivsd
!
!
FUNCTION log_prop(n,x,mean,var_inv)
IMPLICIT none
INTEGER, INTENT(IN) :: n
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x, mean
DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var_inv
DOUBLE PRECISION :: log_prop
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: xm
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: xmtv
ALLOCATE(xm(n),xmtv(n))
xm=x-mean
xmtv=matmul(xm,var_inv)
log_prop=-0.5d0*sum(xmtv*xm)
DEALLOCATE(xm,xmtv)
END FUNCTION log_prop
FUNCTION log_target(nfreq,ncoly,nbeta,x,y,par,sigalpha,tau)
IMPLICIT none
INTEGER, INTENT(IN) :: nfreq,ncoly,nbeta
DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: x
DOUBLE PRECISION, DIMENSION(nfreq,ncoly), INTENT(IN) :: y
DOUBLE PRECISION, DIMENSION(nbeta), INTENT(IN) :: par
DOUBLE PRECISION, INTENT(IN) :: sigalpha, tau
DOUBLE PRECISION :: log_target
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: xpar
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: ys
INTEGER :: j
ALLOCATE(xpar(nfreq),ys(nfreq,ncoly))
xpar=matmul(x,par)
DO j=1,ncoly
ys(:,j)=y(:,j)-xpar
END DO
log_target=sum(ys(2:(nfreq-1),:) - exp(ys(2:(nfreq-1),:)))+ &
0.5d0*sum(ys(1,:)-exp(ys(1,:))) + &
0.5d0*sum(ys(nfreq,:)-exp(ys(nfreq,:)))- &
0.5d0*par(1)**2/sigalpha-0.5d0*sum(par(2:nbeta)**2)/tau
DEALLOCATE(xpar,ys)
END FUNCTION log_target
!
!
FUNCTION log_like(nfreq,nchunks,nbeta,nexp,xx,y,beta,prior_ind)
USE wrrrn_int
IMPLICIT none
INTEGER, INTENT(IN) :: nfreq,nchunks,nbeta,nexp
DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: xx
DOUBLE PRECISION, DIMENSION(nfreq,nchunks), INTENT(IN) :: y
DOUBLE PRECISION, DIMENSION(nexp,nbeta), INTENT(IN) :: beta
DOUBLE PRECISION, DIMENSION(nexp,nchunks), INTENT(IN) :: prior_ind
DOUBLE PRECISION :: log_like
DOUBLE PRECISION, DIMENSION(nfreq,nchunks,nexp) :: ys
DOUBLE PRECISION, DIMENSION(nfreq,nexp) :: xb
DOUBLE PRECISION, DIMENSION(nfreq,nchunks,nexp) :: log_like_arr
DOUBLE PRECISION, DIMENSION(nfreq,nchunks) :: like_mat
INTEGER :: i,j,k
!print*,'log_like func'
!call wrrrn('beta',beta,nra=nexp)
!call wrrrn('xx',xx,nra=nfreq)
DO j=1,nexp
xb(:,j)=matmul(xx,beta(j,:))
END DO
DO j=1,nexp
DO k=1,nchunks
ys(:,k,j)=y(:,k)-xb(:,j)
END DO
END DO
!IF(nexp==1) THEN
! call wrrrn('ys',ys(:,:,1),nra=nfreq)
!END IF
DO j=1,nexp
log_like_arr(2:nfreq-1,:,j)=ys(2:nfreq-1,:,j)-exp(ys(2:nfreq-1,:,j))
log_like_arr(1,:,j)=0.5d0*(ys(1,:,j)-exp(ys(1,:,j)))
log_like_arr(nfreq,:,j)=0.5d0*(ys(nfreq,:,j)-exp(ys(nfreq,:,j)))
END DO
!call wrrrn('log_like_arr',log_like_arr(:,:,1),nra=nfreq)
like_mat=0.0d0
DO k=1,nchunks
DO j=1,nexp
like_mat(:,k)=like_mat(:,k)+prior_ind(j,k)*exp(log_like_arr(:,k,j))
END DO
END DO
log_like=sum(log(like_mat))
END FUNCTION log_like
!
!
FUNCTION log_target_delta(ndelta,nchunks,nexperts,xdelta,z,delta,sigdelta)
USE show_int
IMPLICIT none
INTEGER, INTENT(IN) :: ndelta,nchunks,nexperts
DOUBLE PRECISION, DIMENSION(nchunks,nexperts), INTENT(IN) :: z
DOUBLE PRECISION, DIMENSION(nchunks,2), INTENT(IN) :: xdelta
DOUBLE PRECISION, DIMENSION(ndelta), INTENT(IN) :: delta
DOUBLE PRECISION, INTENT(IN) :: sigdelta
DOUBLE PRECISION :: log_target_delta
DOUBLE PRECISION, DIMENSION(nchunks,nexperts) :: prob
DOUBLE PRECISION, DIMENSION(nchunks,nexperts) :: temp
INTEGER :: i,j,k
temp=0.0d0
DO j=2,nexperts
temp(:,j)=matmul(xdelta,delta(1+(j-2)*2:2*(j-1)))
END DO
DO j=1,nexperts
prob(:,j)=exp(temp(:,j))/sum(exp(temp),2)
END DO
log_target_delta=sum(z*log(prob))-0.5d0*sum(delta**2)/sigdelta
END FUNCTION log_target_delta
!
!
FUNCTION mvrnorm(iseed,n,mean,var)
USE lin_eig_self_int
USE linear_operators
IMPLICIT none
INTEGER, INTENT(IN) :: n,iseed
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: mean
DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var
DOUBLE PRECISION, DIMENSION(n) :: mvrnorm
DOUBLE PRECISION, DIMENSION(n) :: ev
DOUBLE PRECISION, DIMENSION(n) :: z
DOUBLE PRECISION, DIMENSION(n,n) :: evc
DOUBLE PRECISION, DIMENSION(n) :: temp
!find eigen value decomposition of var
call lin_eig_self(var,ev,v=evc)
!draw n independent N(0,1)
call drnnor(n,z)
temp=evc .x. diag(sqrt(ev)) .x. z
mvrnorm=mean+temp
END FUNCTION mvrnorm
!
!
FUNCTION mvrnorm1(iseed,n,mean,var)
USE linear_operators
USE show_int
IMPLICIT none
INTEGER, INTENT(IN) :: n,iseed
DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: mean
DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var
DOUBLE PRECISION, DIMENSION(n) :: mvrnorm1
DOUBLE PRECISION, DIMENSION(n) :: z
DOUBLE PRECISION, DIMENSION(n) :: temp
DOUBLE PRECISION, DIMENSION(n,n) :: fac
DOUBLE PRECISION :: rcond
INTEGER :: i,j
!Find Cholesky decomposition
call dlfcds(n,var,n,fac,n,rcond)
DO i=2,n
DO j=1,i-1
fac(i,j)=0.0d0
END DO
END DO
!draw n independent N(0,1)
call drnnor(n,z)
temp=fac.tx.z
mvrnorm1=mean+temp
END FUNCTION mvrnorm1
!
FUNCTION mvrnorm2(iseed,n,mean,var) RESULT(r)
!n is the number of draws from the multivariate normal
USE rnmvn_int
USE chfac_int
IMPLICIT none
INTEGER, PARAMETER :: DP = kind(1.d0)
INTEGER :: iseed, n
REAL(DP) :: mean(:), var(:,:)
REAL (DP) :: rsig(size(mean),size(mean)), r(n,size(mean))
INTEGER :: irank
INTEGER :: i
call chfac(var,irank,rsig)
call rnmvn(rsig,r)
DO i=1,n
r(i,:)=r(i,:)+mean
END DO
END FUNCTION mvrnorm2
!
!FUNCTION euclid(n,x)
! IMPLICIT none
! INTEGER, INTENT(IN) :: n
! DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x
!END FUNCTION euclid