Posts by author P. Willaume PhiMECA

Module Open TURNS " Integral Compound Poisson Distribution"

Thanks to Anne Dutfoy, a new Open TURNS module is now available. It aims to create a compound Poisson distribution. It is the sum of a Poisson-distributed number of independent identically-distributed random variables.

I invite you to read the documentation shared with the sources of the module to understand the theory and how to use it.

Sources are viewable here: http://trac.openturns.org/browser/openturns-modules/integralcompoundpoisson/trunk

Svn repository is here: https://svn.openturns.org/openturns-modules/integralcompoundpoisson/trunk

Amélioration du chaos - Postprocessing - « Pretty print » des résultats

Ce travail porte sur une présentation lisible des résultats principaux d’un calcul par chaos polynomial : affichage des moments statistiques de la réponse, coefficients du développement en fonction des indices des polynômes, indices de Sobol’.

La fonction printchaos.py (cf. attachments) permet de réaliser un affichage clair des différents résultats du chaos. Elle contient les fonctions :

  • print_moments : dont l'argument est un objet FunctionalChaosResult. En sortie, la fonction affiche la valeur des moments.
import printchaos as pc
...
polynomialChaosAlgorithm.run()
polynomialChaosResult = polynomialChaosAlgorithm.getResult()
pc.print_moments(polynomialChaosResult)
  • print_coefficients : dont l'argument est un objet FunctionalChaosResult. En option, il est possible de spécifier l'objet ComposedDistribution qui a servi à réaliser le calcul par chaos. En sortie, la fonction affiche le type de polynôme relatif à chaque variable si l'objet ComposedDistribution a été spécifié :
    • He : polynôme d'Hermite
    • Le : polynôme de Legendre
    • Ja : polynôme de Jacobi
    • La : polynôme de Laguerre

La fonction affiche ensuite l'ensemble des degrés des polynômes unidimensionnels de chaque variable formant la base de développement du chaos et la valeur du coefficient associée.

import printchaos as pc
...
polynomialChaosAlgorithm.run()
polynomialChaosResult = polynomialChaosAlgorithm.getResult()
pc.print_coefficients(polynomialChaosResult, PM = ProbaModel)
  • print_plot_sensitivities : dont l'argument est un objet FunctionalChaosResult. En option, il est possible de spécifier l'objet ComposedDistribution qui a servi à réaliser le calcul par chaos si les variables aléatoires ont été nommées (méthode setName). La fonction affiche la liste de tous les indices de Sobol’ groupés par ordre en précisant à quelles variables ils sont associés. Enfin la fonction trace un diagramme circulaire des indices de Sobol’ totaux normalisés.
import printchaos as pc
...
polynomialChaosAlgorithm.run()
polynomialChaosResult = polynomialChaosAlgorithm.getResult()
pc.print_plot_sensitivities(polynomialChaosResult, PM = ProbaModel)

Sobol.png, janv. 2011: /extra/oldshare/AmeliorationChaos/Sobol.png

  • print_plot_regression : dont les arguments sont l’échantillon de sortie (NumericalSample) obtenu par à la main (OT 0.13.2) ou par la méthode getProjectionStrategy().getOutputSample() (OT version supérieure) et l’objet FunctionalChaosResult associé au plan d’expérience. La fonction trace un « scatter plot » des résultats du plan d’expérience avec en abscisse les valeurs du vrai modèle et en ordonnées les valeurs du méta-modèle. Dans le terminal, les coefficients de détermination R2 et Q2 sont donnés.
import printchaos as pc
...
polynomialChaosAlgorithm.run()
polynomialChaosResult = polynomialChaosAlgorithm.getResult()
out = polynomialChaosAlgorithm.getProjectionStrategy().getOutputSample()
pc.print_plot_regression(out, polynomialChaosResult)

Model_vs_MetaModel.png, janv. 2011: /extra/oldshare/AmeliorationChaos/Model_vs_MetaModel.png

Amélioration du chaos - Postprocessing - Obtention de la distribution de la réponse d’un chaos polynomial

Il serait utile de pouvoir obtenir, pour la réponse du chaos, un objet Distribution, soit à partir d'un FunctionalChaosResult, soit d'un FunctionalChaosRandomVector. L'objectif est de pouvoir utiliser la distribution de la sortie comme une entrée pour un autre chaos, ou calcul en général. Cette réponse n’est possible qu'avec un lissage à noyaux (KernelSmoothing).

La fonction kernelchaos.py (cf. attachments) permet de créer la distribution de la variable de sortie développée sur le chaos. Elle contient la fonction getOutputDistribution dont l’argument est un objet FunctionalChaosResult. En sortie la fonction retourne une liste de distribution :

  • la distribution lissée par un noyau gaussien,
  • La meilleure distribution usuelle d’après le test de Kolmogorov si le test est validé.
import kernelchaos as kc
...
polynomialChaosAlgorithm.run()
polynomialChaosResult = polynomialChaosAlgorithm.getResult()
ouputDistributions = kc.getOutputDistribution(polynomialChaosResult)

En outre, elle crée 3 graphiques :

  • Output_distribution_kernel_smoothing.png : présentant un échantillon de la variable de sortie sous la forme d’un histogramme ainsi que le lissage à noyau de l’échantillon

Output_distribution_kernel_smoothing.png|C|Output_distribution_kernel_smoothing.png, janv. 2011: /extra/oldshare/AmeliorationChaos/Output_distribution_kernel_smoothing.png

  • Output_distribution_Kolmogorov_test.png : présentant un échantillon de la variable de sortie sous la forme d’un histogramme ainsi que la distribution usuelle validant au mieux le test de Kolmogorov, dans le cas où le test est validé pour un « factory ».

Output_distribution_Kolmogorov_test.png|C|Output_distribution_Kolmogorov_test.png, janv. 2011: /extra/oldshare/AmeliorationChaos/Output_distribution_Kolmogorov_test.png

  • Output_distribution_KK.png : présentant la comparaison entre la distribution lissée par noyau et la distribution usuelle résultat du test de Kolmogorov si ce dernier est validé pour un « factory ».

Output_distribution_KK.png, janv. 2011 /extra/oldshare/AmeliorationChaos/Output_distribution_KK.png

Uncertainity propagation and reliability analysis calculation functions

I made general calculation scripts based on a general model of input data. The syntax of the input data has to follow the following model of script (Attachments: model.py) :

Model of script for input models

  1. The user has to define several model parameters :
  • The number of input random variables;
  • The number of ouput random variables;
  • Finite differente step factor: in Open TURNS 0.12.3 the finite difference step is equal for every variable. In order to take into account the unit of the random variables, I use a finite difference step based on the mean value of each random variable and I define the step as: dx = k_fd*mean(x).

Others parameters can be defined.

...
# --------------------------------------------------------------
# Model parameters
# --------------------------------------------------------------

#Number of input random variables (NbVA)
NbVA=4
#Number of output variables of the physical model (NbVS)
NbVS=1
#Finite differente step factor
k_fd=0.001
#...


#Other parameters
binf=0.0
#...
...
  1. The user has to define his physical model thanks to :
  • wrappers ;
  • OpenTURNSPythonFunction which can use wrappers.

The user has to name his function : PhysModel.

...
# --------------------------------------------------------------
# Physical model
# --------------------------------------------------------------

#Declaration of wrapper(s)
#PhysModel=NumericalMathFunction("my_wrapper")

#Use of OpenTURNSPythonFunction
class myfunction(OpenTURNSPythonFunction):

	def __init__(self):
		OpenTURNSPythonFunction.__init__(self, NbVA, NbVS)

	def f(self, X):
		d=5.0/384.0*X[3]*X[0]**4/(X[1]*X[2])
		return [d]

PhysModel = NumericalMathFunction(myfunction())
...
  1. The user has to define his probabilistic model following several recommandations :
  • each marginal has to be named thanks to the method setName();
  • I only use the independent copula.

The user has to name his composed distribution : ProbaModel.

...
# --------------------------------------------------------------
# Probabilistic model
# --------------------------------------------------------------

# 1. definition of the marginals (1D distribution)

L0=LogNormal(5000.0,50.0,binf,LogNormal.MUSIGMA)
L0.setName('L')
L1=LogNormal(30000.0,4500.0,binf,LogNormal.MUSIGMA)
L1.setName("E")
L2=LogNormal(1.0e9,1.0e8,binf,LogNormal.MUSIGMA)
L2.setName("I")
L3=LogNormal(10.0,3.0,binf,LogNormal.MUSIGMA)
L3.setName("p")
#...

# 2. definition of the distribution collection

CollectionD=DistributionCollection(NbVA)
CollectionD[0]=Distribution(L0)
CollectionD[1]=Distribution(L1)
CollectionD[2]=Distribution(L2)
CollectionD[3]=Distribution(L3)
#...

# 3. Copula definition

IC = IndependentCopula(NbVA)

# 4. Composed distribution definition

ProbaModel = ComposedDistribution(CollectionD,Copula(IC))

# 5. Mean point definition

mean=ProbaModel.getMean()
...
  1. The user can define the types of calculation performed at the beginning of the script (Calculation 1/2). The corresponding functions will run thanks to the second part of the calculation definition (Calculation 2/2). The user has to define the character string related to the numbers of the calculation.
...
# --------------------------------------------------------------
# Calculation 1/2
# --------------------------------------------------------------
# 0: Single evaluation of the physical model at mean point

# 11: Quadratic cumul analysis
# 12: Uncertainity propagation using Monte Carlo simulations

# 121: Basic statistics on Monte Carlo results

# in the future
# 13: Polynomial chaos expansion

# 21: Reliability analysis using FORM
# 22: Reliability analysis using SORM

# 31: Reliability analysis using Monte Carlo simulations
# 32: Reliability analysis using latin hypercube sampling
# 33: Reliability analysis using directional simulations

# 211: Reliability analysis using importance sampling


Calculation=["0","12","121"]

for i in range(len(Calculation)):
	if (Calculation[i]=="11"):
		from fc_CQ import Calcul_CQ
	if (Calculation[i]=="12"):
		from fc_MC import Calcul_MC

	if (Calculation[i]=="121"):
		from fc_BasicStat import Calcul_BasicStat

...

# --------------------------------------------------------------
# Calculation 2/2
# --------------------------------------------------------------

#Modification of the finite difference step value
step=NumericalPoint(NbVA)
for i in range(NbVA):
	step[i]=mean[i]*k_fd
myGradient = NonCenteredFiniteDifferenceGradient(step,PhysModel.getEvaluationImplementation())
myHessian = CenteredFiniteDifferenceHessian(step,PhysModel.getEvaluationImplementation())
PhysModel.setGradientImplementation(myGradient)
PhysModel.setHessianImplementation(myHessian)

for i in range(len(Calculation)):
	if (Calculation[i]=="0"):
		Mean_evaluation=PhysModel(mean)
		for j in range(NbVS):
			print ""
			print "Evaluation of the physical model at mean point - output",j+1,"/",NbVS,":", Mean_evaluation[j]
	if (Calculation[i]=="11"):
		resCQ=Calcul_CQ(PhysModel,ProbaModel)
	if (Calculation[i]=="12"):
		resMC=Calcul_MC(PhysModel,ProbaModel)

	if (Calculation[i]=="121"):
		resBasicStat=Calcul_BasicStat(resMC)
...

The calculation functions are independent files, except Calcul_BasicStat. Each file contains a python function based on Open TURNS. Each function can be divided into 2 parts :

  • the calculation definition and run;
  • the display of the results.

Function Quadratic cumul analysis

Attachements: fc_CQ.py

The function requires:

  • a NumericalMathFunction or physical model
  • a ComposedDistribution or probabilistic model

Note : In Open TURNS 0.12.3, the calculation of the importance factors is only available for scalar NumericalMathFunction.

...
# Quadratic cumul analysis
	RV_ProbabilisticModel = RandomVector(Distribution(ProbabilisticModel))
	Function = RandomVector(PhysicalModel,RV_ProbabilisticModel)
	QuadraticC = QuadraticCumul(Function)

	if (NbVS==1):
		IF=NumericalPoint(NbVA)
		for i in range(NbVA):
			IF[i] = QuadraticC.getImportanceFactors()[i]

	Mean1=QuadraticC.getMeanFirstOrder()
	Mean2=QuadraticC.getMeanSecondOrder()

	CovG=QuadraticC.getCovariance()

	stdev=NumericalPoint(NbVS)
	for i in range(NbVS):
		stdev[i] = sqrt(CovG[i,i])
...

Function Uncertainity propagation using Monte Carlo simulations

Attachements: fc_MC.py

The function requires:

  • a NumericalMathFunction or physical model
  • a ComposedDistribution or probabilistic model
...
# Monte Carlo
	In = ProbabilisticModel.getNumericalSample(10000)
	Out = PhysicalModel(In)

	Moyenne = Out.computeMean()
	standarddeviation=Out.computeStandardDeviationPerComponent()
	Skewness = Out.computeSkewnessPerComponent()
	Kurtosis = Out.computeKurtosisPerComponent()
...

Function Basic statistics

Attachements: fc_BasicStat.py

The function requires:

  • a NumericalSample (result of uncertainity propagation using Monte Carlo simulations.

This kind of calculation has to be performed just after an uncertainity propagation using Monte Carlo simulations. This function will be able to be improved.

Last modifications : 12/06/2009

First script

you can test Open TURNS with this simple script :

from openturns import *
 
dist = Normal( 0, 1 )
sample = dist.getNumericalSample(10)
print sample

Now you have printed a sample !

Enjoy.