1#! /usr/bin/env python 2 3from __future__ import print_function 4import openturns as ot 5 6ot.TESTPREAMBLE() 7 8# Instanciate one distribution object 9dim = 3 10R = ot.CorrelationMatrix(dim) 11for i in range(dim - 1): 12 R[i, i + 1] = 0.25 13copula = ot.SklarCopula( 14 ot.Distribution(ot.Normal([1.0, 2.0, 3.0], [2.0, 3.0, 1.0], R))) 15copulaRef = ot.NormalCopula(R) 16print("Copula ", repr(copula)) 17print("Copula ", copula) 18print("Mean =", repr(copula.getMean())) 19print("Mean (ref)=", repr(copulaRef.getMean())) 20ot.ResourceMap.SetAsUnsignedInteger("GaussKronrod-MaximumSubIntervals", 20) 21ot.ResourceMap.SetAsScalar("GaussKronrod-MaximumError", 1.0e-4) 22print("Covariance =", repr(copula.getCovariance())) 23ot.ResourceMap.SetAsUnsignedInteger("GaussKronrod-MaximumSubIntervals", 100) 24ot.ResourceMap.SetAsScalar("GaussKronrod-MaximumError", 1.0e-12) 25print("Covariance (ref)=", repr(copulaRef.getCovariance())) 26 27# Is this copula an elliptical distribution? 28print("Elliptical distribution= ", copula.isElliptical()) 29 30# Is this copula elliptical ? 31print("Elliptical copula= ", copula.hasEllipticalCopula()) 32 33# Is this copula independent ? 34print("Independent copula= ", copula.hasIndependentCopula()) 35 36# Test for realization of distribution 37oneRealization = copula.getRealization() 38print("oneRealization=", repr(oneRealization)) 39 40# Test for sampling 41size = 10 42oneSample = copula.getSample(size) 43print("oneSample=", repr(oneSample)) 44 45# Test for sampling 46size = 1000 47anotherSample = copula.getSample(size) 48print("anotherSample mean=", repr(anotherSample.computeMean())) 49print("anotherSample covariance=", repr(anotherSample.computeCovariance())) 50 51# Define a point 52point = ot.Point(dim, 0.2) 53 54# Show PDF and CDF of point 55pointPDF = copula.computePDF(point) 56pointCDF = copula.computeCDF(point) 57pointPDFRef = copulaRef.computePDF(point) 58pointCDFRef = copulaRef.computeCDF(point) 59print("Point = ", repr(point), " pdf =%.6f" % 60 pointPDF, " cdf =%.6f" % pointCDF) 61print("Point = ", repr(point), " pdf (ref)=%.6f" % 62 pointPDFRef, " cdf (ref)=%.6f" % pointCDFRef) 63 64# Get 50% quantile 65quantile = copula.computeQuantile(0.5) 66quantileRef = copulaRef.computeQuantile(0.5) 67print("Quantile =", repr(quantile)) 68print("Quantile (ref)=", repr(quantileRef)) 69print("CDF(quantile)=%.6f" % copula.computeCDF(quantile)) 70# Get 95% survival function 71inverseSurvival = ot.Point(copula.computeInverseSurvivalFunction(0.95)) 72print("InverseSurvival=", repr(inverseSurvival)) 73print("Survival(inverseSurvival)=%.6f" % 74 copula.computeSurvivalFunction(inverseSurvival)) 75# Takes too much time 76# print("entropy=%.6f" % copula.computeEntropy()) 77# Confidence regions 78if copula.getDimension() <= 2: 79 threshold = ot.Point() 80 print("Minimum volume interval=", 81 copula.computeMinimumVolumeInterval(0.95, threshold)) 82 print("threshold=", threshold) 83 beta = ot.Point() 84 levelSet = copula.computeMinimumVolumeLevelSet(0.95, beta) 85 print("Minimum volume level set=", levelSet) 86 print("beta=", beta) 87 print("Bilateral confidence interval=", 88 copula.computeBilateralConfidenceInterval(0.95, beta)) 89 print("beta=", beta) 90 print("Unilateral confidence interval (lower tail)=", 91 copula.computeUnilateralConfidenceInterval(0.95, False, beta)) 92 print("beta=", beta) 93 print("Unilateral confidence interval (upper tail)=", 94 copula.computeUnilateralConfidenceInterval(0.95, True, beta)) 95 print("beta=", beta) 96 97# Extract the marginals 98for i in range(dim): 99 margin = copula.getMarginal(i) 100 marginRef = copulaRef.getMarginal(i) 101 print("margin=", repr(margin)) 102 print("margin PDF =%.6f" % 103 margin.computePDF(ot.Point(1, 0.25))) 104 print("margin PDF (ref)=%.6f" % 105 marginRef.computePDF(ot.Point(1, 0.25))) 106 print("margin CDF =%.6f" % 107 margin.computeCDF(ot.Point(1, 0.25))) 108 print("margin CDF (ref)=%.6f" % 109 marginRef.computeCDF(ot.Point(1, 0.25))) 110 print("margin quantile =", repr(margin.computeQuantile(0.95))) 111 print("margin quantile (ref)=", repr(marginRef.computeQuantile(0.95))) 112 print("margin realization=", repr(margin.getRealization())) 113 114# Extract a 2-D marginal 115indices = ot.Indices([1, 0]) 116print("indices=", repr(indices)) 117margin = copula.getMarginal(indices) 118marginRef = copulaRef.getMarginal(indices) 119print("margin=", repr(margin)) 120print("margin PDF =%.6f" % margin.computePDF(ot.Point(2, 0.25))) 121print("margin PDF (ref)=%.6f" % 122 marginRef.computePDF(ot.Point(2, 0.25))) 123print("margin CDF =%.6f" % margin.computeCDF(ot.Point(2, 0.25))) 124print("margin CDF (ref)=%.6f" % 125 marginRef.computeCDF(ot.Point(2, 0.25))) 126print("margin quantile =", repr(margin.computeQuantile(0.95))) 127print("margin quantile (ref)=", repr(marginRef.computeQuantile(0.95))) 128print("margin realization=", repr(margin.getRealization())) 129 130# tbb nested parallelism issue 131student = ot.Student(3.0, [1.0]*2, [3.0]*2, ot.CorrelationMatrix(2)) 132copula = ot.SklarCopula(student) 133pdf_graph = copula.drawPDF() 134cdf_graph = copula.drawCDF() 135