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