Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

Register Log in
  • Congratulations strongm on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!

Multivariate normal pdf code?

Status
Not open for further replies.
ignacio82

ignacio82

Technical User
Feb 15, 2012
9
US
I'm looking for Fortran code that calculates the multivariate normal pdf.

If someone is willing to share or point me in the right direction I would appreciate it.

Thanks!
 
Sort by date Sort by votes
If you haven't found anything on google and it is an engineering type application, you could try the sister site www.eng-tips.com. They are more clued up about where to find stuff like that.
 
Upvote 0 Downvote
sorry this reply is so late

please find below some fortran code that possibly does what you want -- multivariate normal pdf.

i worked on a small project several years ago. part of it involved calculating the multivariate normal pdf. this is a sample of my code. it should help you get started. sorry the code is a bit messy / sloppy :(

please let me know if you can use it

i hope it helps


PS sorry i copied / pasted the code, instead of uploading it (although there's not much code). i know it's not the etiquette. but, i had a little difficulty with uploading. i might have better luck next time :)



Code: 
      SUBROUTINE CALC_XPDF_X1X2X3(X,IX,NX,MEAN_X,IMEAN_X,SD_X,ISD_X,COV_X,&
						&ICOV_X,JCOV_X,INV_COV_X,IINV_COV_X,JINV_COV_X,&
						&DET_COV_X,RHO_X12,RHO_X13,RHO_X23,XPDF,INTEST)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!		  
!     PURPOSE:
!     --------
!
!     CALC_XPDF_X1X2X3 computes the (joint PDF) of X.
!
!
!     X - Vector of variables (Normal distribution) where fX(x) is to be
!         - determined
!     XPDF - Output as fX(x)
!     MEAN_X,SD_X - Mean and Standard Deviation of X, respectively
!                 - ie Marginal PDF for ith variable (Xi)
!     COV_X - Normal Covariance Matrix for X
!     INV_COV_X - Inverse of COV_X
!     DET_COV_X - Determinant of COV_X
!     RHO_XIJ - Correlation coefficient of XI,XJ
!     INTEST   -  Error checking flag
!             -  Contains a value corresponding to the following errors
!             -  0   No errors
!             -  1   XPDF <0.0 
!             -  2   XPDF can't be computed
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      INTEGER,INTENT(IN) ::IX,NX,IMEAN_X,ISD_X,ICOV_X,JCOV_X
      INTEGER,INTENT(IN) ::IINV_COV_X,JINV_COV_X
	  REAL(KIND=HIGH),INTENT(IN),DIMENSION(IX) :: X
	  REAL(KIND=HIGH),INTENT(IN),DIMENSION(IMEAN_X) :: MEAN_X
	  REAL(KIND=HIGH),INTENT(IN),DIMENSION(ISD_X) :: SD_X
	  REAL(KIND=HIGH),INTENT(IN),DIMENSION(ICOV_X,JCOV_X) :: COV_X
	  REAL(KIND=HIGH),INTENT(IN),DIMENSION(IINV_COV_X,JINV_COV_X) :: INV_COV_X
	  REAL(KIND=HIGH),INTENT(IN) :: DET_COV_X
	  REAL(KIND=HIGH),INTENT(IN) :: RHO_X12,RHO_X13,RHO_X23
	  REAL(KIND=HIGH),INTENT(OUT) :: XPDF
      INTEGER,INTENT(INOUT) ::INTEST


      INTEGER ::IERR,I
	  REAL(KIND=HIGH) :: PDF
	  REAL(KIND=HIGH),PARAMETER :: SQRT_TWO_PI=2.506628275D0
    

      IF (NX==1)	THEN

		CALL PDF_NORM(X(1),MEAN_X(1),SD_X(1),XPDF,INTEST)

      ELSE IF ((RHO_X12==0.0D0).AND.(RHO_X13==0.0D0).AND.(RHO_X23==0.0D0))	THEN

		XPDF=1.0D0
		DO I=1,TOTLOAD
			CALL PDF_NORM(X(I),MEAN_X(I),SD_X(I),PDF,INTEST)
			XPDF=XPDF*PDF
		END DO

      ELSE 

		CALL JOINT_PDF_NORM_NX3(X,IX,MEAN_X,IMEAN_X,SD_X,ISD_X,COV_X,&
							&ICOV_X,JCOV_X,INV_COV_X,IINV_COV_X,&
							&JINV_COV_X,DET_COV_X,NX,XPDF,INTEST)

      END IF





      END SUBROUTINE CALC_XPDF_X1X2X3
!********************************************************************
      SUBROUTINE PDF_NORM(X,XM,XSD,PDFX,INTEST)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!		  
!     PURPOSE:
!     --------
!
!     PDF_NORM computes the PDF at X, for a NORMal distribution.
!
!
!     X - Location where PDF is to be computed
!     XM - Mean of X
!     XSD - Standard deviation of X
!     PDFX - Output as fx(X)
!     INTEST   -  Error checking flag
!             -  Contains a value corresponding to the following errors
!             -  0   No errors
!             -  1   PDFX <0.0 
!
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	  REAL(KIND=HIGH),INTENT(IN) :: X,XM,XSD
	  REAL(KIND=HIGH),INTENT(OUT) :: PDFX
      INTEGER,INTENT(INOUT) ::INTEST


      INTEGER ::IERR
	  REAL(KIND=HIGH) :: U
	  REAL(KIND=HIGH),PARAMETER :: SQRT_TWO_PI=2.506628275D0
    

	  U=(X-XM)/XSD
	  PDFX=DEXP(-0.5D0*U*U)
	  PDFX=PDFX/SQRT_TWO_PI/XSD




      END SUBROUTINE PDF_NORM
!********************************************************************
      SUBROUTINE JOINT_PDF_NORM_NX3(X,IX,XM,IXM,XSD,IXSD,COV_X,&
							&ICOV_X,JCOV_X,INV_COV_X,IINV_COV_X,&
							&JINV_COV_X,DET_COV_X,NX,PDFX,INTEST)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!		  
!     PURPOSE:
!     --------
!
!     JOINT_PDF_NORM_NX3 computes the JOINT PDF at X, for a NORMal
!     distribution, for a system having >=3 variables (NX>=3).
!
!
!     X - Vector Location where PDF is to be computed
!     XM - Vector of Mean of X
!     XSD - Vector of Standard Deviation of X
!     COV_X - Normal Covariance Matrix for X
!     INV_COV_X - Inverse of COV_X
!     DET_COV_X - Determinant of COV_X
!     NX - Number of random variables
!     PDFX - Output as joint fx(X)
!     INTEST   -  Error checking flag
!             -  Contains a value corresponding to the following errors
!             -  0   No errors
!             -  1   PDFX <0.0 
!
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      INTEGER,INTENT(IN) ::IX,NX,IXM,IXSD,ICOV_X,JCOV_X
      INTEGER,INTENT(IN) ::IINV_COV_X,JINV_COV_X
	  REAL(KIND=HIGH),INTENT(IN),DIMENSION(IX) :: X
	  REAL(KIND=HIGH),INTENT(IN),DIMENSION(IXM) :: XM
	  REAL(KIND=HIGH),INTENT(IN),DIMENSION(IXSD) :: XSD
	  REAL(KIND=HIGH),INTENT(IN),DIMENSION(ICOV_X,JCOV_X) :: COV_X
	  REAL(KIND=HIGH),INTENT(IN),DIMENSION(IINV_COV_X,JINV_COV_X) :: INV_COV_X
	  REAL(KIND=HIGH),INTENT(IN) :: DET_COV_X
	  REAL(KIND=HIGH),INTENT(OUT) :: PDFX
      INTEGER,INTENT(INOUT) ::INTEST


      INTEGER ::IERR
      INTEGER ::I,J,N,M
	  REAL(KIND=HIGH) :: X1,X2,SD_X1,SD_X2,MEAN_X1,MEAN_X2,RHO,H,K
	  REAL(KIND=HIGH) :: U
	  REAL(KIND=HIGH) :: SUM_AB,SUM_ABC
	  REAL(KIND=HIGH),DIMENSION(IX) :: X_XM
	  REAL(KIND=HIGH),DIMENSION(IX) :: AB
	  REAL(KIND=HIGH),DIMENSION(IX) :: X_XMT
	  REAL(KIND=HIGH),PARAMETER :: TWO_PI=6.283185307D0
    


!     compute multivariate PDF

!     compute (X-XM) and (X-XM)T

		  X_XM=0.0D0
		  X_XMT=0.0D0
		  DO I=1,NX
			X_XM(I)=X(I)-XM(I)
			X_XMT(I)=X(I)-XM(I)
		  END DO

!     compute (X-XM)T*INVERSE(COV_Q)
	  

			AB=0.0D0
			DO J=1,NX
				SUM_AB=0.0D0
				DO N=1,NX
					SUM_AB=SUM_AB+X_XMT(N)*INV_COV_X(N,J)
				END DO
				AB(J)=SUM_AB
			END DO


!     compute (X-XM)T*INVERSE(COV_X)*(X-XM)

		  SUM_ABC=0.0D0
		  DO M=1,NX
			  SUM_ABC=SUM_ABC+AB(M)*X_XM(M)
		  END DO

		   U=SUM_ABC

!     compute PDFX

		   PDFX=DEXP(-0.5D0*U)
		   PDFX=PDFX*DET_COV_X**(-0.5)
		   PDFX=PDFX*TWO_PI**(-DBLE(NX)/2.0D0)





      END SUBROUTINE JOINT_PDF_NORM_NX3
 
Upvote 0 Downvote
Status
Not open for further replies.

Similar threads

Space_Fede
  • Locked
  • Question
Fortran code to read unformatted binary files 1
2
Replies
34
Views
574
JamesONeil
JamesONeil
ASRC8367
  • Locked
  • Question
Fortran code to average hourly data into daily values 1
Replies
2
Views
132
mikrom
mikrom
I_rvln
  • Locked
  • Question
Launch fortran code on linux to read data
Replies
6
Views
160
mikrom
mikrom
Scarletttt
  • Locked
  • Question
Help with compiling Fortran Legacy code
Replies
1
Views
106
xwb
xwb
FloatingFeather
  • Locked
  • Question
Fortran90 code with FORTRAN77 subroutine 1
Replies
3
Views
88
xwb
xwb

Part and Inventory Search

Sponsor

Back
Top