Blog: Amélioration du chaos - Postprocessing - Obtention de la distribution de la réponse d’un chaos polynomial: kernelchaos.py

File kernelchaos.py, 3.6 KB (added by souchaud@…, 8 years ago)
Line 
1import openturns as ot
2
3def getOutputDistribution(res):
4
5    N1 = 1000
6    N2 = 50
7
8    myRandomVector_O = ot.RandomVector(ot.FunctionalChaosRandomVector(res))
9
10    myNumericalSample_O1 = myRandomVector_O.getNumericalSample(N1)
11    kernelNormal = ot.Distribution(ot.Normal(0.0,1.0))
12    kernel_O = ot.KernelSmoothing(ot.Distribution(kernelNormal))
13    myBandwidth_O = kernel_O.computeSilvermanBandwidth(myNumericalSample_O1)
14
15    print ""
16    print "Generation of the output distribution by kernel smoothing"
17    print ""
18
19    myKernelSmoothing_O = kernel_O.buildImplementation(myNumericalSample_O1, myBandwidth_O)
20    myDistribution_O1 = ot.Distribution(myKernelSmoothing_O)
21
22
23    plot_O1 = myDistribution_O1.drawPDF()
24
25    hist=ot.VisualTest.DrawHistogram(myNumericalSample_O1)
26    plot_hist=hist.getDrawables()[0]
27
28    plot_O1.addDrawable(ot.Drawable(plot_hist))
29    plot_O1.setTitle("Gaussian kernel smoothing on\n"+str(N1)+" simulations of the meta model")
30    plot_O1.setLegendPosition("")
31    plot_O1.draw("Output_distribution_kernel_smoothing", 800, 600, ot.GraphImplementation.PNG)
32
33
34
35    myDistributionFactoryCollection = ot.DistributionFactoryCollection(9)
36    myDistributionFactoryCollection[0] = ot.DistributionFactory(ot.BetaFactory())
37    myDistributionFactoryCollection[1] = ot.DistributionFactory(ot.ExponentialFactory())
38    myDistributionFactoryCollection[2] = ot.DistributionFactory(ot.GammaFactory())
39    myDistributionFactoryCollection[3] = ot.DistributionFactory(ot.GumbelFactory())
40    myDistributionFactoryCollection[4] = ot.DistributionFactory(ot.LogisticFactory())
41    myDistributionFactoryCollection[5] = ot.DistributionFactory(ot.LogNormalFactory())
42    myDistributionFactoryCollection[6] = ot.DistributionFactory(ot.NormalFactory())
43    myDistributionFactoryCollection[7] = ot.DistributionFactory(ot.RayleighFactory())
44    myDistributionFactoryCollection[8] = ot.DistributionFactory(ot.WeibullFactory())
45
46
47    myNumericalSample_O2 = myRandomVector_O.getNumericalSample(N2)
48    myDistribution_O2 = ot.FittingTest.BestModelKolmogorov(myNumericalSample_O2, myDistributionFactoryCollection)
49    myFittingTestResult_O2 = ot.FittingTest.Kolmogorov(myNumericalSample_O2, myDistribution_O2)
50
51
52    print 'Kolmogorov p-value =', myFittingTestResult_O2.getPValue()
53    print 'Test succeed ?', myFittingTestResult_O2.getBinaryQualityMeasure()
54    print myDistribution_O2.getName(), ":", myDistribution_O2.getParametersCollection()[0]
55
56
57    if myFittingTestResult_O2.getBinaryQualityMeasure():
58        plot_O2 = myDistribution_O2.drawPDF()
59
60        hist=ot.VisualTest.DrawHistogram(myNumericalSample_O2)
61        plot_hist=hist.getDrawables()[0]
62
63        plot_O2.addDrawable(ot.Drawable(plot_hist))
64        plot_O2.setTitle("Kolmogorov fitting test on\n"+str(N2)+" simulations of the meta model")
65        plot_O2.setLegendPosition("")
66        plot_O2.draw("Output_distribution_Kolmogorov_test", 800, 600, ot.GraphImplementation.PNG)
67
68        plot_O1O2 = myDistribution_O1.drawPDF()
69        plot=plot_O1O2.getDrawables()[0]
70        plot.setLegendName("Gaussian kernel smoothing")
71
72        O22 = myDistribution_O2.drawPDF()
73        plot_O22=O22.getDrawables()[0]
74        plot_O22.setColor("blue")
75        plot_O22.setLegendName("Kolmogorov distribution:\n"+myDistribution_O2.getName()+":"+str(myDistribution_O2.getParametersCollection()[0]))
76
77        plot_O1O2.addDrawable(ot.Drawable(plot_O22))
78        plot_O1O2.setTitle("Kernel vs. Kolmogorov")
79        plot_O1O2.draw("Output_distribution_KK", 800, 600, ot.GraphImplementation.PNG)
80
81
82
83        return [myDistribution_O1, myDistribution_O2]
84    else:
85        return [myDistribution_O1]
86