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