MODULE CODE4
  USE CODE5, ONLY : Sample_Mixture_Normal,Sample_Gamma,Sample_truncated_Normal, &
       Sample_Multivariate_Normal,Sample_Double_Truncated_Normal,NORMCDF,MIXNORMPDF, &
       Sample_Uniform
  USE CODE9
  USE CODE6
  USE CODE7, ONLY : MATINV
  IMPLICIT NONE

CONTAINS

  REAL(8) FUNCTION SAMPLE_F(this,mf,tf)     
    IMPLICIT NONE
    TYPE(UNIVARIATE_NORMAL_MIXTURE), INTENT(IN) :: this
    REAL(8), INTENT(IN) :: mf,tf
    REAL(8) :: p(this%k),t(this%k),m(this%k),aux(this%k)
    INTEGER :: k
    DO k = 1, this%k
       t(k) = tf + this%tau(k)
       m(k) = (tf*mf + this%tau(k)*this%mu(k))/t(k)
       aux(k) = -0.5d0*(tf*this%tau(k)/t(k))*((this%mu(k)-mf)**2)
       p(k) = this%l(k)*DSQRT(tf*this%tau(k)/t(k))*DEXP(aux(k))
    END DO
    p = p / SUM(p)
    SAMPLE_F = Sample_Mixture_Normal(m,1.0d0/t,p)    
  END FUNCTION SAMPLE_F

  ! samples tau 
  REAL(8) FUNCTION TAU(g1,g2,Y,X,f,b,a)
    REAL(8), INTENT(IN) :: g1,g2,Y(:),X(:,:),b(:),f(:,:),a(:)
    REAL(8) :: a1,a2
    a1 = 0.5d0*DBLE(SIZE(Y,1)) + g1
    a2 = 0.5d0*SUM((Y-MATMUL(X,b)-MATMUL(f,a))**2) + g2
    TAU = Sample_gamma(a1,a2)
  END FUNCTION TAU

  SUBROUTINE SAMPLE_ALPHA(this,prior_var)
    IMPLICIT NONE
    REAL(8), INTENT(IN) :: prior_var
    TYPE(linear_factor_equation), INTENT(INOUT) :: this
    REAL(8) :: P(this%nor%k_f,this%nor%k_f),a(this%nor%k_f)
    INTEGER :: j
    P = this%t*MATMUL(TRANSPOSE(this%f_f),this%f_f)
    DO j = 1, this%nor%k_f
       P(j,j) = P(j,j) + prior_var
    END DO
    P=MATINV(P)
    a = this%t*MATMUL(TRANSPOSE(this%f_f),this%yq)
    a = MATMUL(P,a)
    this%a(this%nor%i_f)=Sample_Multivariate_Normal(a,P)
  END SUBROUTINE SAMPLE_ALPHA

  SUBROUTINE SAMPLE_ALPHA_TRUNCATED(this,prior_var)
    IMPLICIT NONE
    REAL(8), INTENT(IN) :: prior_var
    TYPE(linear_factor_equation), INTENT(INOUT) :: this
    REAL(8) :: P,a
    INTEGER :: j
    P = this%t*DOT_PRODUCT(this%f_s,this%f_s)
    P = P + prior_var
    P = 1.0d0/P
    a = this%t*DOT_PRODUCT(this%f_s,this%yq)
    a = P*a
    this%a(this%nor%i_s) = Sample_Truncated_Normal(a,P,0.0d0,.true.)
  END SUBROUTINE SAMPLE_ALPHA_TRUNCATED

  SUBROUTINE SAMPLE_BETA(this)
    IMPLICIT NONE
    TYPE(linear_factor_equation), INTENT(INOUT) :: this
    REAL(8) :: P(this%k,this%k),b(this%k)
    INTEGER :: j
    P = this%ixx/this%t
    b = MATMUL(TRANSPOSE(this%x),this%yq)
    b = MATMUL(this%ixx,b)
    this%b = Sample_Multivariate_Normal(b,P)
  END SUBROUTINE SAMPLE_BETA

END MODULE CODE4
