1#! /usr/bin/env python
2
3from __future__ import print_function, division
4import unittest as ut
5from math import pi, log, log10, floor
6
7import openturns as ot
8# (value, reference, rtol, atol)
9from openturns.testing import assert_almost_equal
10
11ot.TESTPREAMBLE()
12
13# Prevent from using unittest.TestLoader.loadTestsFromTestCase
14# so as not to get the test execution duration
15
16
17def loadTestsFromTestCase(cls):
18    """Launch the tests of the test case class"""
19    import sys
20    from inspect import getmembers
21    members = getmembers(cls)
22    test_names = []
23    is_tearDownClass_there = False
24    for member in members:
25        member_name = member[0]
26        if member_name[0:5] == 'test_':
27            test_names.append(member_name)
28        elif member_name == 'setUpClass':
29            cls.setUpClass()
30        elif member_name == 'tearDownClass':
31            is_tearDownClass_there = True
32    for test_name in test_names:
33        test = cls(test_name)
34        test.setUp()
35        print("Run " + test_name + "... ", end="")
36        sys.stdout.flush()
37        try:
38            test.debug()
39            print("SUCCESS")
40        except Exception as exception:
41            print("FAILURE")
42            print(exception)
43        test.tearDown()
44    if is_tearDownClass_there:
45        cls.tearDownClass()
46
47
48class TestInverseWishartMethods(ut.TestCase):
49
50    """Test case for the class InverseWishart"""
51
52    @classmethod
53    def setUpClass(cls):
54        # attributes to compare a one-dimensional InverseWishart to the
55        # equivalent InverseGamma distribution
56        U = ot.Uniform(0., 1.)
57        scale = 10. * U.getRealization()[0]
58        DoF = 3. + U.getRealization()[0]  # Degrees of Freedom
59        cls.k, cls.beta = 0.5 * DoF, 0.5 * scale
60        cls.one_dimensional_inverse_wishart \
61            = ot.InverseWishart(ot.CovarianceMatrix([[scale]]), DoF)
62        cls.inverse_gamma = ot.InverseGamma(1. / cls.beta, cls.k)
63        # attributes to test a multi-dimensional InverseWishart
64        cls.dimension = 5
65        cls.DoF = cls.dimension + 3 + U.getRealization()[0]
66        cls.L = ot.TriangularMatrix(cls.dimension)
67        diagL = ot.Uniform(0.5, 1.).getSample(cls.dimension)
68        cls.determinant = 1.
69        for i in range(cls.dimension):
70            cls.determinant *= diagL[i, 0]
71            cls.L[i, i] = diagL[i, 0]
72            for j in range(i):
73                cls.L[i, j] = U.getRealization()[0]
74        cls.determinant *= cls.determinant
75        cls.Scale = ot.CovarianceMatrix(cls.L * cls.L.transpose())
76
77    def test_computeLogPDF_1D_case(self):
78        """Test InverseWishart.computeLogPDF in the one-dimensional case"""
79        k, beta = self.k, self.beta
80
81        def logPDF(x):
82            if x <= 0.:
83                raise ValueError("math domain error")
84            return k * log(beta) - ot.SpecFunc_LogGamma(k) - (k + 1) * log(x) - beta / x
85        data = ((self.inverse_gamma.drawPDF()).getDrawable(0)).getData()
86        i = 0
87        while data[i, 0] <= 0.:
88            i += 1
89        for d in data[i:, 0]:
90            x = d[0]
91            logPDFx = logPDF(x)
92            logPDFx_IW = self.one_dimensional_inverse_wishart.computeLogPDF(x)
93            logPDFx_IG = self.inverse_gamma.computeLogPDF(x)
94            assert_almost_equal(logPDFx_IW, logPDFx)
95            assert_almost_equal(logPDFx_IG, logPDFx)
96            assert_almost_equal(logPDFx_IW, logPDFx_IG)
97
98    # Not a test
99    # The log multi-gamma function appears in the log PDF
100    def logmultigamma(self, p, x):
101        """Computes the logarithm of the multi-gamma function at x"""
102        logmgpx = 0.25 * p * (p - 1) * log(pi)
103        for i in range(1, p + 1):
104            logmgpx = logmgpx + ot.SpecFunc_LnGamma(x + 0.5 * (1 - i))
105        return logmgpx
106
107    # Test InverseWishart.computeLogPDF in the special case of diagonal matrices
108    # (scale covariance matrix and point/matrice at which to evaluate the
109    # PDF) by comparing the logarithm of the ratio of the InverseWishart PDF
110    # by the product of the InverseGamma PDFs
111    def test_computeLogPDF_diagonal_case(self):
112        """Test InverseWishart.computeLogPDF in the case of diagonal matrices"""
113        dimension, DoF = self.dimension, self.DoF
114        k = 0.5 * (DoF + dimension - 1)
115        diagX = ot.Uniform(0.5, 1.).getSample(dimension)
116        Scale = ot.CovarianceMatrix(dimension)
117        X = ot.CovarianceMatrix(dimension)
118        for d in range(dimension):
119            Scale[d, d], X[d, d] = self.Scale[d, d], diagX[d, 0]
120        inverse_wishart = ot.InverseWishart(Scale, DoF)
121        logdensity = inverse_wishart.computeLogPDF(X)
122        logratio = - self.logmultigamma(dimension, 0.5 * DoF) \
123            + dimension * ot.SpecFunc_LnGamma(0.5 * (DoF + dimension - 1))
124        for d in range(dimension):
125            inverse_gamma = ot.InverseGamma(2. / Scale[d, d], k)
126            logdensity = logdensity - inverse_gamma.computeLogPDF(diagX[d, 0])
127            logratio = logratio + 0.5 * \
128                (1 - dimension) * log(0.5 * Scale[d, d])
129        assert_almost_equal(logdensity, logratio)
130
131    # Test InverseWishart.computeLogPDF by evaluating the log PDF
132    # at the scale covariance matrix
133    def test_computeLogPDF(self):
134        """Test InverseWishart.computeLogPDF"""
135        dimension, DoF = self.dimension, self.DoF
136        Scale, determinant = self.Scale, self.determinant
137        inverse_wishart = ot.InverseWishart(Scale, DoF)
138        logPDFatX = - self.logmultigamma(dimension, 0.5 * DoF) \
139            - 0.5 * (DoF * dimension * log(2.) + dimension
140                     + (dimension + 1) * log(determinant))
141        assert_almost_equal(inverse_wishart.computeLogPDF(Scale), logPDFatX)
142
143    # Compare the empirical expectations of a large matrix sample
144    # and of the corresponding inverse matrix sample
145    # to the corresponding theoretical expectations
146    def test_getSample_getMean(self):
147        """Test InverseWishart.getSample and InverseWishart.getMean"""
148        d, Scale, DoF, N = self.dimension, self.Scale, self.DoF, int(1E+4)
149        Identity = ot.CovarianceMatrix(d)
150        Scale_wishart = ot.CovarianceMatrix(Scale.solveLinearSystem(Identity))
151        inverse_wishart = ot.InverseWishart(Scale, DoF)
152        sample_inverse = ot.Sample(N, (d * (d + 1)) // 2)
153        sample = ot.Sample(N, (d * (d + 1)) // 2)
154        for i in range(N):
155            M_inverse = inverse_wishart.getRealizationAsMatrix()
156            M = M_inverse.solveLinearSystem(Identity)
157            indice = 0
158            for j in range(d):
159                for k in range(j + 1):
160                    sample_inverse[i, indice] = M_inverse[k, j]
161                    sample[i, indice] = M[k, j]
162                    indice += 1
163        mean_inverse = sample_inverse.computeMean()
164        mean = sample.computeMean()
165        theoretical_mean_inverse = inverse_wishart.getMean()
166        theoretical_mean = (ot.Wishart(Scale_wishart, DoF)).getMean()
167        indice, coefficient = 0, 1. / (DoF - d - 1)
168        for j in range(d):
169            for k in range(j + 1):
170                assert_almost_equal(theoretical_mean_inverse[indice],
171                                    coefficient * Scale[k, j])
172                assert_almost_equal(theoretical_mean[indice],
173                                    DoF * Scale_wishart[k, j])
174                assert_almost_equal(mean_inverse[indice],
175                                    coefficient * Scale[k, j], 0.15, 1.E-3)
176                assert_almost_equal(mean[indice],
177                                    DoF * Scale_wishart[k, j], 0.15, 1.E-3)
178                indice += 1
179
180
181# Instanciate one distribution object
182distribution = ot.InverseWishart(ot.CovarianceMatrix(1), 5.0)
183print("Distribution ", repr(distribution))
184print("Distribution ", distribution)
185
186# Get mean and covariance
187print("Mean= ", repr(distribution.getMean()))
188print("Covariance= ", repr(distribution.getCovariance()))
189
190# Is this distribution elliptical ?
191print("Elliptical = ", distribution.isElliptical())
192
193# Test for realization of distribution
194oneRealization = distribution.getRealization()
195print("oneRealization=", repr(oneRealization))
196
197# Test for sampling
198size = 10000
199oneSample = distribution.getSample(size)
200print("oneSample first=", repr(
201    oneSample[0]), " last=", repr(oneSample[size - 1]))
202print("mean=", repr(oneSample.computeMean()))
203print("covariance=", repr(oneSample.computeCovariance()))
204
205size = 100
206for i in range(2):
207    msg = ''
208    if ot.FittingTest.Kolmogorov(distribution.getSample(size), distribution).getBinaryQualityMeasure():
209        msg = "accepted"
210    else:
211        msg = "rejected"
212    print("Kolmogorov test for the generator, sample size=", size, " is", msg)
213    size *= 10
214
215# Define a point
216point = ot.Point(distribution.getDimension(), 9.1)
217print("Point= ", repr(point))
218
219# Show PDF and CDF of point
220eps = 1e-5
221# max = distribution.getB() + distribution.getA()
222# min = distribution.getB() - distribution.getA()
223# derivative of PDF with regards its arguments
224DDF = distribution.computeDDF(point)
225print("ddf     =", repr(DDF))
226# by the finite difference technique
227print("ddf (FD)=", repr(ot.Point(1, (distribution.computePDF(point + ot.Point(1, eps)) -
228                                     distribution.computePDF(point + ot.Point(1, -eps))) / (2.0 * eps))))
229
230# PDF value
231LPDF = distribution.computeLogPDF(point)
232print("log pdf=%.6f" % LPDF)
233PDF = distribution.computePDF(point)
234print("pdf     =%.6f" % PDF)
235# by the finite difference technique from CDF
236print("pdf (FD)=%.6f" % ((distribution.computeCDF(point + ot.Point(1, eps)) -
237                          distribution.computeCDF(point + ot.Point(1, -eps))) / (2.0 * eps)))
238
239# derivative of the PDF with regards the parameters of the distribution
240CDF = distribution.computeCDF(point)
241print("cdf=%.6f" % CDF)
242CCDF = distribution.computeComplementaryCDF(point)
243print("ccdf=%.6f" % CCDF)
244PDFgr = distribution.computePDFGradient(point)
245print("pdf gradient     =", repr(PDFgr))
246# by the finite difference technique
247PDFgrFD = ot.Point(2)
248v00 = distribution.getV()[0, 0]
249nu = distribution.getNu()
250PDFgrFD[0] = (ot.InverseWishart(ot.CorrelationMatrix([[v00]]), nu).computePDF(point) -
251              ot.InverseWishart(ot.CorrelationMatrix([[v00 - eps]]), nu).computePDF(point)) / (1.0 * eps)
252PDFgrFD[1] = (ot.InverseWishart(ot.CorrelationMatrix([[v00]]), nu + eps).computePDF(point) -
253              ot.InverseWishart(ot.CorrelationMatrix([[v00]]), nu - eps).computePDF(point)) / (2.0 * eps)
254print("pdf gradient (FD)=", repr(PDFgrFD))
255
256# derivative of the PDF with regards the parameters of the distribution
257CDFgr = distribution.computeCDFGradient(point)
258print("cdf gradient     =", repr(CDFgr))
259CDFgrFD = ot.Point(2)
260CDFgrFD[0] = (ot.InverseWishart(ot.CorrelationMatrix([[v00]]), nu).computeCDF(point) -
261              ot.InverseWishart(ot.CorrelationMatrix([[v00 - eps]]), nu).computeCDF(point)) / (1.0 * eps)
262CDFgrFD[1] = (ot.InverseWishart(ot.CorrelationMatrix([[v00]]), nu + eps).computeCDF(point) -
263              ot.InverseWishart(ot.CorrelationMatrix([[v00]]), nu - eps).computeCDF(point)) / (2.0 * eps)
264print("cdf gradient (FD)=",  repr(CDFgrFD))
265
266# quantile
267quantile = distribution.computeQuantile(0.95)
268print("quantile=", repr(quantile))
269print("cdf(quantile)=%.6f" % distribution.computeCDF(quantile))
270# Get 95% survival function
271inverseSurvival = ot.Point(distribution.computeInverseSurvivalFunction(0.95))
272print("InverseSurvival=", repr(inverseSurvival))
273print("Survival(inverseSurvival)=%.6f" %
274      distribution.computeSurvivalFunction(inverseSurvival))
275print("entropy=%.6f" % distribution.computeEntropy())
276
277# Confidence regions
278interval, threshold = distribution.computeMinimumVolumeIntervalWithMarginalProbability(
279    0.95)
280print("Minimum volume interval=", interval)
281print("threshold=", ot.Point(1, threshold))
282levelSet, beta = distribution.computeMinimumVolumeLevelSetWithThreshold(0.95)
283print("Minimum volume level set=", levelSet)
284print("beta=", ot.Point(1, beta))
285interval, beta = distribution.computeBilateralConfidenceIntervalWithMarginalProbability(
286    0.95)
287print("Bilateral confidence interval=", interval)
288print("beta=", ot.Point(1, beta))
289interval, beta = distribution.computeUnilateralConfidenceIntervalWithMarginalProbability(
290    0.95, False)
291print("Unilateral confidence interval (lower tail)=", interval)
292print("beta=", ot.Point(1, beta))
293interval, beta = distribution.computeUnilateralConfidenceIntervalWithMarginalProbability(
294    0.95, True)
295print("Unilateral confidence interval (upper tail)=", interval)
296print("beta=", ot.Point(1, beta))
297
298mean = distribution.getMean()
299print("mean=", repr(mean))
300standardDeviation = distribution.getStandardDeviation()
301print("standard deviation=", repr(standardDeviation))
302skewness = distribution.getSkewness()
303print("skewness=", repr(skewness))
304kurtosis = distribution.getKurtosis()
305print("kurtosis=", repr(kurtosis))
306covariance = distribution.getCovariance()
307print("covariance=", repr(covariance))
308parameters = distribution.getParametersCollection()
309print("parameters=", repr(parameters))
310for i in range(6):
311    print("standard moment n=", i, " value=",
312          distribution.getStandardMoment(i))
313print("Standard representative=", distribution.getStandardRepresentative())
314
315loadTestsFromTestCase(TestInverseWishartMethods)
316