Does imsl for fortran have a multivariate normal cdf routine?
I found the instructions for c but not for fortran.
Thanks!
Does imsl for fortran have a multivariate normal cdf routine?
I found the instructions for c but not for fortran.
Thanks!
Fortran does not have a multivariate normal cdf routine. But C does, as you note. If you have access to the C library, the Fortran 2003 standard supports an intrinsic module for inter-operating. This means the Fortran code that you write can call the C routine. I have used this feature several times in my own research. Writing this interface code is a short investment in time and certainly beats waiting for the code to appear in a Fortran library. So the key question for you is whether you have a license to use the C library.
Last edited by Richard Hanson; 10-06-2012 at 11:06 AM.
Dickie T. Bird
I have a licence to use the C library. Could you point me the right direction about how to do this?
Thanks!
I will work up an example and post it later. Other readers may benefit, so it is worth while to go to some work. What you will see has these elements:
- You define a Fortran derived type that packages the optional arguments that the C routine can accept.
- Match the default value settings of the C routine by initialization statements within the derived type.
- Give the derived type a C binding when it is defined. This fixes the order of the components and other details.
- Provide a Fortran interface to a C user-wrapper function that will unpack the derived type argument (now considered a C structure in the wrapper code).
- The Fortran code then calls the wrapper code, using the C binding. The derived type is included as an argument.
- This C wrapper code is in a small library, created as a separate item.
- The wrapper code calls the multivariate CDF code with all optional arguments and returns the value to Fortran.
- The linking of the Fortran code will include the IMSL C statistic library and the small additional library.
There are ways to short-cut these steps, and the style is personal. But this outline can be used to add additional C routines for use by the Fortran programmer. So stand by. Don't worry if this is confusing. The example should clear things up. After that it is up to you for getting your compilation and linking steps correct!
Last edited by Richard Hanson; 10-08-2012 at 10:58 AM.
Dickie T. Bird
The three files attached illustrate the inter-language call of Fortran calling the IMSL C stat code for the multi-variate normal CDF. Understand that this is not supported by IMSL. Please do not ask for support. It is provided as an illustration of how to use the Fortran 2003 standard and call C routines from Fortran. There is a lot more to this subject. For example see the white paper on the Rogue Wave website about this topic.
Standardized Mixed Language Programming for Fortran and C
http://www.roguewave.com/documents.aspx?EntryId=797
Last edited by Richard Hanson; 10-08-2012 at 12:01 PM.
Dickie T. Bird
Thanks a lot!
I am coding in fortran and using your code above to calculate the multivariate cdf by calling the C library. I got the code to work fine. I got the correct answer for your example as well as the examples on the imsl website. However, when I started entering in cdf's that I needed to calculate weird things started happening. In particular, I want to find the cdf of a normal variable with
mean= [0,.75,0]
variance covariance matrix=
[1,0,0]
[0,1,-1]
[0,-1,2].
The upper bounds for calculating the cdf that I am using are:
[1.5,1.5,0].
When I do this using your code, I get cdf=1.0. When I do this in matlab I get .2791
If I change the variance-covariance matrix to:
[1,0,0]
[0,1,-.9999]
[0,-.9999,1]
Then I get that the cdf=.2551 using your code; matlab gets that the cdf=.2791
Lastly, if I make the variance covariance matrix
[1,0,0]
[0,1,-1.0001]
[0,-1.0001,1]
I get the cdf=1 using your code and equals .2791 using matlab.
I attached my code below.
I would really appreciate any advice you could offer me in why I am getting such drastic changes in the cdf when the variance-covariance matrix is slightly changed.
Thanks!
Last edited by cookicat53; 10-15-2012 at 10:38 AM.
I am violating my admonition to the poster: This posted code is an illustration only and shows how to call the C routine from Fortran. No support is intended. If you feel the C routine is incorrect then contact post-sales support staff. But your issue is an inconsistency in the results. My pair of matched compilers (Intel/MS C++) does not show any such problem. See the output below. In your posting you did not have the most SE element at 2. always, but the code does.
From looking at this output, the issue now is that Matlab gives you the result .279 and the IMSL code gives you .255, both approximately. Explaining this difference has to go back to the algorithm. IMSL uses a method by Alan Genz. Perhaps Matlab does, too. Clearly there are some differences. What they are will be hard to determine unless both source codes are available to scan or else a correct result is available from some other source. Sorry for your trouble. This sort of thing is really frustrating. So get in touch with post-sales support if and when you find wrong results in the IMSL C code.
Code:Upper limits of integration: 1.5000 1.5000 0.0000 Means in each dimension: 0.0000 0.7500 0.0000 Correlation matrix: 1.0000 0.0000 0.0000 0.0000 1.0000 -1.0000 0.0000 -1.0000 2.0000 Probability and error estimate: 0.255109 1.0000E-12 Upper limits of integration: 1.5000 1.5000 0.0000 Means in each dimension: 0.0000 0.7500 0.0000 Correlation matrix: 1.0000 0.0000 0.0000 0.0000 1.0000 -0.9999 0.0000 -0.9999 2.0000 Probability and error estimate: 0.255109 1.0000E-12 Upper limits of integration: 1.5000 1.5000 0.0000 Means in each dimension: 0.0000 0.7500 0.0000 Correlation matrix: 1.0000 0.0000 0.0000 0.0000 1.0000 -1.0001 0.0000 -1.0001 2.0000 Probability and error estimate: 0.255109 1.0000E-12
Last edited by Richard Hanson; 10-15-2012 at 05:51 PM.
Dickie T. Bird
This thread has remained unattended to for six years, with an apparent discrepancy between the results from IMSL and from Matlab (or other alternatives). I happened to revisit this thread this week, and I am pleased to explain the reason for the discrepancy and to provide a simple solution to resolve that discrepancy.
I wish to pause to note that Dr. Richard J. Hanson, who spent many years working on IMSL until he left Roguewave in 2013, died of cancer in December 2016; an obituary written by his friend, collaborator and well-known mathematician Dr. Fred Krogh (author of the Mathalacarte web site and the JPL MATH77 library) is available at https://sinews.siam.org/Details-Page...chard-j-hanson . Hanson was one of the creators of BLAS. Rest in Peace, Dr. Hanson!
The CNL routines imsls_f_multivariate_normal_cdf and imsls_d_multivariate_normal_cdf take as their 4th argument a two-dimensional array, sigma. The CNL 8.0 documentation ambiguously describes this argument as "Array of length k by k containing the positive definite symmetric matrix of correlations or of variances and covariances for the multivariate normal distribution". This needs to be corrected to say "Array of length k by k containing the positive definite symmetric matrix of correlations". The relation between the correlation matrix and the variance/covariance matrix is available at, for example, https://en.wikipedia.org/wiki/Covari...nd_correlation . The (i,j) element of the covariance matrix should be divided by the square root of the product of the (i,i) and the (j,j) elements of the same matrix in order to obtain the (i,j) element of the correlation matrix.
For the example problem described by Cookicat53 in #7, the input lines for the correlation matrix should be:
1.0, 0.0, 0.0,
0.0, 1.0, -0.707106781186548,
0.0, -0.707106781186548, 1.0
When Hanson's Fortran 90 code is corrected this way, CNL 8.0 returns the correct answer, which is 0.279074 .This result agrees with the Matlab result quoted by Cookicat53, and with the result given by the NAG FL library routine G01HBF. That routine effectively re-scales the input sigma-matrix so that the main diagonal contains all elements = 1, and the CNL routine could do the same internally. If the CNL 8.0/8.6 routine were so modified, the CNL documentation could be left as it is.
Last edited by mecej4; 01-02-2019 at 09:59 AM.