Chaospy: An open source tool for designing methods of uncertainty quantification

ContentslistsavailableatScienceDirect

Journal of Computational Science

jo u r n al ho me p a g e :w w w . e l s e v i e r . c o m / l o c a t e / j o c s

Chaospy: An open source tool for designing methods of uncertainty quantification

Jonathan Feinberg

, Hans Petter Langtangen

aCenterforBiomedicalComputing,SimulaResearchLaboratory,P.O.Box134,Lysaker,Norway

bDepartmentofMathematics,UniversityofOslo,P.O.Box1053,Blindern,Oslo,Norway

cDepartmentofInformatics,UniversityofOslo,P.O.Box1080,Blindern,Oslo,Norway

a r t i c l e i n f o

Articlehistory:

Received19March2015

Receivedinrevisedform19July2015 Accepted16August2015

Availableonline29August2015

Keywords:

Uncertaintyquantification Polynomialchaosexpansions MonteCarlosimulation Rosenblatttransformations Pythonpackage

a b s t r a c t

Thepaperdescribesthephilosophy,design,functionality,andusageofthePythonsoftwaretoolbox ChaospyforperforminguncertaintyquantificationviapolynomialchaosexpansionsandMonteCarlo simulation.ThepapercomparesChaospytosimilarpackagesanddemonstratesastrongerfocusondefin- ingreusablesoftwarebuildingblocksthatcaneasilybeassembledtoconstructnew,tailoredalgorithms foruncertaintyquantification.Forexample,aChaospyusercaninafewlinesofhigh-levelcomputercode definecustomdistributions,polynomials,integrationrules,samplingschemes,andstatisticalmetricsfor uncertaintyanalysis.Inaddition,thesoftwareintroducessomenovelmethodologicaladvances,likea frameworkforcomputingRosenblatttransformationsandanewapproachforcreatingpolynomialchaos expansionswithdependentstochasticvariables.

©2015TheAuthors.PublishedbyElsevierB.V.ThisisanopenaccessarticleundertheCCBYlicense (http://creativecommons.org/licenses/by/4.0/).

1. Introduction

Weconsideracomputationalscienceprobleminspace xand timetwheretheaimistoquantifytheuncertaintyinsomeresponse Y,computedbyaforwardmodelf,whichdependsonuncertain inputparameters Q:

Y=f(x,t,Q). (1)

Wetreat Q asavectorofmodelparameters,andYisnormally computedassomegridfunctioninspaceandtime.Theuncertainty inthisproblemstemsfromtheparameters Q,whichareassumed tohaveaknownjointprobabilitydensityfunctionp_Q.Thechallenge isthatwewanttoquantifytheuncertainty inY,butnothingis knownaboutitsdensityp_Y.Thegoalisthentoeitherbuildthe densitypYorrelevantdescriptivepropertiesofYusingthedensity p_Qandtheforwardmodelf.Forallpracticalpurposesthismustbe donebyanumericalprocedure.

Inthispaper,wefocusontwoapproachestonumericallyquan- tifyuncertainty:MonteCarlosimulationandnon-intrusiveglobal polynomialchaosexpansions.Forareviewoftheformer,thereis averyusefulbookbyRubinstein,ReuvenandKroese[1],whilefor thelatter,werefertotheexcellentbookbyXiu[2].Notethatother methodsforperforminguncertaintyquantificationalsoexist,such

∗ Correspondingauthorat:CenterforBiomedicalComputing,SimulaResearch Laboratory,P.O.Box134,Lysaker,Norway.

asperturbationmethods,momentequations,andoperatorbased methods.Thesemethodsarealldiscussedin[2],butarelessgeneral andlesswidelyapplicablethanthetwoaddressedinthispaper.

The number of toolboxes available to perform Monte Carlo simulationisvastlylargerthanthenumberoftoolboxesfornon- intrusivepolynomialchaosexpansion.Asfarastheauthorsknow, thereareonlyafewviableoptionsforthelatterclassofmethods:

TheDakotaProject(referredtoasDakota)[3],theOpusOpenTurns library(referredtoasTurns)[4],UncertentyQuantificationToolkit [5],andMITUncertentyQuantificationLibrary[6].Inthispaperwe willfocusontheformer two:Dakota andTurns.Bothpackages consistoflibrarieswithextensivesetsoftools,whereMonteCarlo simulationandnon-intrusivepolynomialchaosexpansionarejust twotoolsavailableamongseveralothers.

ItisworthnotingthatbothDakotaandTurnscanbeusedfrom twoperspectives:asauserandasadeveloper.Bothpackagesare opensourceprojectswithcomprehensivedevelopermanuals.As such,theybothallowanyonetoextendthesoftwarewithanyfunc- tionalityoneseesfit.However,theseextensionfeaturesarenot targetingthecommonuserandrequireadeeperunderstandingof bothcodingpracticeandtheunderlyingdesignofthelibrary.In ouropinion,thethresholdforacommonusertoextendthelibrary isnormallyoutofreach.Consequently,weareinthispaperonly consideringDakotaandTurnsfromthepointofviewofthecommon user.

Dakotarequirestheforwardmodelftobewrappedinastand- alonecallableexecutable.Thecommonapproach isthentolink

http://dx.doi.org/10.1016/j.jocs.2015.08.008

1877-7503/©2015TheAuthors.PublishedbyElsevierB.V.ThisisanopenaccessarticleundertheCCBYlicense(http://creativecommons.org/licenses/by/4.0/).

thisexecutabletotheanalysissoftwarethroughaconfiguration file.Thetechnicalstepsaresomewhatcumbersome,buthastheir advantageinthatalreadybuiltandinstalledsimulationsoftware canbeusedwithoutwritingalineofcode.

Alternativetothisdirectapproachistointeractwithanapplica- tionprogramminginterface(API).Thisapproachrequirestheuser toknowhowtoprograminthesupportedlanguages,butthisalso hasclearbenefitsasaninterfacethroughaprogramminglanguage allowsforadeeperlevelofintegrationbetweentheuser’smodel andtheUQtools.Also,exposingthesoftware’sinternalcomponents throughanAPIallowsahigherdetailedcontroloverthetoolsand howtheycanbecombinedinstatisticalalgorithms.Thisfeatureis attractivetoscientistswhowouldlikethepossibilitytoexperiment withnewornon-standardmethodsinwaysnotthoughtofbefore.

ThisapproachisusedbytheTurnssoftware(usingthelanguages PythonorR)andissupportedinDakotathroughalibrarymode (usingC++).

For example, consider bootstrapping [7], a popular method formeasuringthestabilityofanyparameterestimation.Neither DakotanorTurnssupportbootstrappingdirectly.However,since Turnsexposessomeoftheinnercomponentstotheuser,apro- grammercancombinethesetoimplementacustombootstrapping technique.

Thispaperdescribesanew,thirdalternativeopensourcesoft- warepackagecalledChaospy [8].LikeDakota andTurns, itis a toolboxforanalysinguncertaintyusingadvancedMonteCarlosim- ulationandnon-intrusivepolynomialchaosexpansions.However, unliketheothers,itaimstoassistscientistsinconstructingtail- oredstatisticalmethodsbycombiningalot offundamentaland advancedbuildingblocks.Chaospybuildsuponthesamephilos- ophyas Turnsin thatit offersflexibility totheuser,but takes itsignificantlyfurther.InChaospy,itispossibletogaindetailed controlandadduserdefinedfunctionalitytoallofthefollowing:

random variablegeneration, polynomial construction,sampling schemes, numerical integration rules, response evaluation, and pointcollocation.Thesoftwareisdesignedfromthegroundupin Pythontobemodularandeasytoexperimentwith.Thenumber oflinesofcodetoachieveafulluncertaintyanalysisisamazingly low.Itisalsoveryeasytocomparearangeofmethodsinagiven problem.Standardstatisticalmethodsareeasilyaccessiblethrough afewlinesofRorPandas[9]code,andonemaythinkofChaospy asatoolsimilartoRorPandas,justtailoredtopolynomialchaos expansionandMonteCarlosimulation.

AlthoughChaospyisdesignedwithalargefocusonmodularity, flexibility,andcustomization,thetoolboxcomeswithawiderange ofpre-definedstatisticalmethods.WithinthescopeofMonteCarlo samplingandnon-intrusivepolynomialchaosexpansion,Chaospy hasacompetitivecollectionofmethods,comparabletobothDakota andTurns.Italsoofferssomenovelfeaturesregardingstatistical methods,firstand foremostlya flexible frameworkfor defining andhandlinginputdistributions,includingdependentstochastic variables.Detailedcomparisonsoffeaturesinthethreepackages appearthroughoutthepaper.

Thepaperisstructuredasfollows.WestartinSection2witha quickdemonstrationofhowthesoftwarecanbeusedtoperform uncertaintyquantificationina simplephysicalproblem.Section 3addressesprobability distributionsandthetheoryrelevantto performMonteCarlosimulation.Section4concernsnon-intrusive polynomialchaosexpansions,whileconclusionsandtopicsforfur- therworkappearinSection5.

2. AglimpseofChaospyinaction

TodemonstratehowChaospyisusedtosolveanuncertainty quantificationproblem,weconsiderasimplephysicalexampleof

(scaled)exponentialdecaywithanuncertain,piecewiseconstant coefficient:

u(x)=−c(x)u(x), u(0)=u0, c(x)=

⎧ ⎪

⎨

⎪ ⎩

c0, x<0.5 c1, 0.5≤x<0.7 c2, x≥0.7

(2)

Suchamodelarisesinmanycontexts,butwemayherethinkofu(x) astheporosityatdepthxingeologicallayersandc_iasa(scaled) compactionconstantinlayernumberi.Forsimplicity,weconsider onlythreelayerswiththreeuncertainconstantsc0,c1,andc2.

Themodelcaneasilybeevaluatedbysolvingthedifferential equationproblem,herebya2nd-orderRunge–Kuttamethodona meshx,codedinPythonas:

Alternatively,themodelcanbeimplementedinsomeexter- nalsoftwareinanotherprogramminglanguage.Thissoftwarecan eitherberunasastand-aloneapplication,wherethePythonfunc- tionmodelrunstheapplicationandcommunicateswithitthrough inputandoutputfiles,orthemodelfunctioncancommunicatewith theexternalsoftwarethroughfunctioncallsifaPythonwrapper hasbeenmadeforthesoftware(there arenumeroustechnolo- giesavailableforcreatingPythonwrappersforC,C++,andFortran software).

TheChaospypackagemaybeloadedby

Eachoftheuncertainparametersmustbeassignedaprobabil- itydensity,andwe assumethatc0,c1,andc2 arestochastically independent:

Thesamplepoints(c₀,c₁,c₂)inprobabilityspace,wherethe modelistobeevaluated,canbechoseninmanyways.Herewe specifyathird-orderGaussianQuadratureschemetailoredtothe jointdistribution:

Thenextstepistoevaluatethecomputationalmodelatthese samplepoints(objectnodes):

Now,samplescontainsalistofarrays,eacharraycontainingu valuesatthe101xvaluesforonecombination(c₀,c₁,c₂)ofthe inputparameters.

Fig. 1. Solution of a simple stochastic differential equation with uncertain coefficients.

To create a polynomial chaos expansion,we must generate orthogonalpolynomialscorrespondingtothejointdistribution.We choosepolynomialsofthesameorderasspecifiedinthequadrature rule,computedbythewidelyusedthree-termrecurrencerelation (ttr):

To create an approximate solver (or surrogate model), we jointhepolynomialchaosexpansion,thequadraturenodesand weights,andthemodelsamples:

Themodelapproxobjectcannowcheaplyevaluatethemodel atapoint(c₀,c₁,c₂)inprobabilityspaceforallxpointsinthex array.Built-intoolscanbeusedtoderivestatisticalinformation aboutthemodelresponse:

Themeananddeviationobjectsarearrayscontainingthemean valueandstandarddeviationateachpointinx.Agraphicalillus- trationisshowninFig.1.

TheaccuracyoftheestimationiscomparabletowhatDakota andTurnscanprovide.Fig.2showsthattheestimationerrorinthe threesoftwaretoolboxesarealmostindistinguishable.Theerroris calculatedastheabsolutedifferencebetweenthetruevalueand theestimatedvalueintegratedoverthedepthx:

εE=

|E(u)−E(uapprox)|dx εV=

|V(u)−V(uapprox)|dx

Boththepointcollocationmethodandthepseudo-spectralprojec- tionmethodareincluded.Theformeriscalculatedusingtwotimes therandomcollocationnodesasthenumberofpolynomials,and thelatterusingGaussianquadratureintegrationwithquadrature orderequaltopolynomialorder.NotethatTurnsdoesnotsupport pseudo-spectralprojection,andisthereforeonlycomparedusing pointcollocation.

3. Modellingrandomvariables 3.1. Rosenblatttransformation

Numericalmethodsforuncertaintyquantificationneedtogen- eratepseudo-randomrealizations

{Qk}_k∈IK IK={1,...,K},

fromthedensityp_Q.Each Q∈{Qk}_k∈IK is multivariatewiththe numberofdimensionsD>1.Generatingrealizationsfromagiven densityp_Qisoftennon-trivial,atleastwhenDislarge.Averycom- monassumptionmadeinuncertaintyquantificationisthateach dimensionin Q consistsof stochasticallyindependentcompo- nents.Stochasticindependenceallowsforajointsamplingscheme tobereducedtoaseriesofunivariatesamplings,drasticallyreduc- ingthecomplexityofgeneratingasample Q.

Unfortunately,theassumptionofindependencedoesnotalways holdin practice.We have examples frommany researchfields wherestochasticdependencemustbeassumed, includingmod- ellingofclimate[10],iron-oreminerals[11],finance[12],andion channeldensitiesindetailedneurosciencemodels[13].Therealso existsexampleswhereintroducingdependentrandomvariablesis beneficialforthemodellingprocess,eventhoughtheoriginalinput wasstochasticallyindependent[14].In anycases, modellingof stochasticallydependentvariablesarerequiredtoperformuncer- taintyquantificationadequately.AstrongfeatureofChaospyisits supportforstochasticdependence.

Allrandom samples arein Chaospy generated using Rosen- blatttransformationsT_Q [15].Itallowsforarandomvariable U,

Fig.2.Theerrorinestimatesofthemeanandvariance,computedbyDakota,Turns,andChaospyusingpointcollocationandpseudo-spectralprojection,isalmostidentical.

generateduniformlyonaunithypercube[0,1]^D,tobetransformed intoQ=T_Q⁻¹(U),whichbehavesasifitweredrawnfromthedensity p_Q.Itiseasytogeneratepseudo-randomsamplesfromauniform distribution,andtheRosenblatttransformationcanthenbeused asamethodforgeneratingsamplesfromarbitrarydensities.

TheRosenblatttransformationcanbederivedasfollows.Con- sider a probability decomposition, for example for a bivariate randomvariable Q=(Q₀,Q₁):

pQ0,Q1(q0,q1)=pQ0(q0)pQ1|Q0(q1|q0), (3) werepQ0isanmarginaldensityfunction,andpQ1|Q0isaconditional density.Forthemultivariatecase,thedensitydecompositionwill havetheform

p_Q(q)=

D−1

d=0

p_Q

d(q_d), (4)

where

Q_d =Qd|Q0,...,Qd−1 q_d=qd|q0,...,qd−1 (5) denotesthatQ_dandQ_daredependentonallcomponentswithlower indices.AforwardRosenblatttransformationcanthenbedefined as

TQ(q)=(FQ₀(q₀),...,FQ_{D− 1}(q_D−1)), (6) whereF_Q

disthecumulativedistributionfunction:

FQ_d(q_d)=

−∞

pQ_d(r|q0,...,qd−1)dr. (7)

Thistransformationisbijective,soitisalwayspossibletodefine theinverseRosenblatttransformationT_Q⁻¹inasimilarfashion.

3.2. NumericalestimationofinverseRosenblatttransformations Toimplement theRosenblatt transformation in practice,we needtoidentifytheinversetransformT_Q⁻¹.Unfortunately,T_Q is oftennon-linearwithouta closed-formformula,makinganalyt- icalcalculations of thetransformation’s inversedifficult.In the scenariowherewedonothaveasymbolicrepresentationofthe inversetransformation,anumericalschemehastobeemployed.

Totheauthors’knowledge,therearenostandardsfordefiningsuch anumericalscheme.Thefollowingparagraphsthereforedescribe ourproposedmethodforcalculating theinversetransformation numerically.

TheproblemofcalculatingtheinversetransformationT_Q⁻¹can, bydecomposingthedefinitionoftheforwardRosenblatttransfor- mationin(6),bereformulatedas

F_Q⁻¹ d

(u|q0,...,qd−1)

=

r:FQ_d(r|q0,...,q_d−1)=u d=0,...,D−1.

Inotherwords,thechallengeofcalculatingtheinversetransforma- tioncanbereformulatedasaseriesofonedimensionalroot-finding problems. In Chaospy, these roots are found by employing a Newton–Raphsonscheme.However,toensureconvergence,the schemeiscoupledwithabisectionmethod.Thebisectionmethod is applicable here since the problem is one-dimensional and thefunctionsof interestareby definitionmonotone.Whenthe Newton–Raphson method fails to converge at an increment, a bisectionstepgivestheNewton–Raphsonanewstartlocationaway fromthepreviouslocation.Thisalgorithmensuresfastandreliable convergencetowardstheroot.

TheNewton–Raphson-bisectionhybridmethodisimplemented asfollows.Theinitialvaluesarethelowerandupperbounds[lo0,

up₀].IfpQ_d isunbound,theintervalisselectedsuchthatitapprox- imately covers thedensity. For example for a standard normal random variable,which isunbound,theinterval [−7.5,7.5] will approximatelycoverthewholedensitywithanerrorabout10⁻¹⁴. ThealgorithmstartswithaNewton–Raphsonincrement,usingthe initialvaluer0=(up0−lo0)u+lo0:

r_k+1=rk−FQ_d(rk|q0,...,qd−1)−u p_Q

d(r_k|q0,...,q_d−1) , (8)

wherethedensitypQ

dcanbeapproximatedusingfinitedifferences.

Ifthenewvaluedoesnotfallintheinterval[lo_k,up_k],thispro- posedvalueisrejected, andisinsteadreplacedwithabisection increment:

r_k+1= upk+lok

2 . (9)

Ineithercase,theboundsareupdatedaccordingto

(lo_k+1,up_k+1)=

d(r_k+1|q0,...,q_d−1)<u (10) The algorithm repeats thesteps in (8)–(10), until theresidual

|F_Q

d(rk|q0,...,q_d−1)−u|issufficientlysmall.

Thedescribed algorithm overcomesoneof thechallengesof implementingRosenblatttransformationsinpractice:howtocal- culatethe inversetransformation. Anotherchallenge is how to constructatransformationinthefirstplace.Thisisthetopicof thenextsection.

3.3. Constructingdistributions

ThebackboneofdistributionsinChaospyistheRosenblatttrans- formationT_Q.Themethod,asdescribedintheprevioussection, assumesthatp_Q isknowntobeabletoperformthetransforma- tionanditsinverse.Inpractice,however,wefirstneedtoconstruct p_Q,beforethetransformation can beused. Thiscanbe a chal- lenging task, but in Chaospy a lot of effort has been putinto constructingnoveltools for makingtheprocess asflexible and painlessaspossible.Inessence,userscancreatetheirowncustom multivariatedistributionsusinganewmethodologyasdescribed next.

Followingthedefinitionin(6),eachRosenblatttransformation consistsofacollectionofconditionaldistributions.Weexpressall conditionalitythroughdistributionparameters.Forexample,the locationparameterofanormaldistributioncanbesettobeuni- formlydistributed,sayon[−1,1].ThefollowinginteractivePython codedefinesanormalvariablewithanormallydistributedmean:

Wenowhavetwostochasticvariables,uniformandnormal, whosejointbivariatedistributioncanbeconstructedthroughthe cp.Jfunction:

Thesoftwarewill,fromthisminimalformulation,trytosortout thedependencyorderingandconstructthefullRosenblatttrans- formation.Theonlyrequirementisthatadecompositionasin(4) is infactpossible.Theresult isa fullyfunctioningforwardand inverseRosenblatttransformation. Thefollowing codeevaluates theforwardtransformation(the density)at(1,0.9),theinverse

Table1

Listofsupportedcontinuousdistributionsinsoftware.Thetitles‘D’,‘T’and‘C’representsDakota,TurnsandChaospyrespectively.Theelements‘y’and‘n’representthe answers‘yes’and‘no’indicatingifthedistributionissupportedornot.

Distribution D T C Distribution D T C

Alpha n n y Anglit n n y

Arcsinus n n y Beta y y y

Brandford n n y Burr n y y

Cauchy n n y Chi n y y

Chi-square n y y DoubleGamma n n y

DoubleWeibull n n y Epanechnikov n y y

Erlang n n y Exponential y y y

Exponentialpower n n y ExponentialWeibull n n y

Birnbaum–Sanders n n y Fisher–Snedecor n y y

Fisk/log-logistic n n y FoldedCauchy n n y

Foldednormal n n y Frechet y n y

Gamma y y y Gen.exponential n n y

Gen.extremevalue n n y Gen.gamma n n y

Gen.half-logistic n n y Gilbrat n n y

TruncatedGumbel n n y Gumbel y y y

Hypergeometricsecant n n y Inverse-normal n y n

Kumaraswamy n n y Laplace n y y

Levy n n y Log-gamma n n y

Log-laplace n n y Log-normal y y y

Log-uniform y y y Logistic n y y

Lomax n n y Maxwell n n y

Mielke’sbeta-kappa n n y Nakagami n n y

Non-centralchi-squared n y y Non-centralStudent-T n y y

Non-centralF n n y Normal y y y

Pareto(firstkind) n n y Powerlog-normal n n y

Powernormal n n y Raisedcosine n n y

Rayleigh n y y Reciprocal n n y

Rice n y n Right-skewedGumbel n n y

Student-T n y y Trapezoidal n y n

Triangle y y y Truncatedexponential n n y

Truncatednormal n y y Tukey-Lambda n n y

Uniform y y y Wald n n y

Weibull y y y Wigner n n y

WrappedCauchy n n y Zipf–Mandelbrot n y n

transformationat(0.4,0.6),anddrawsarandomsamplefromthe jointdistribution:

Distributionsin higher dimensions are triviallyobtained by includingmoreargumentstothecp.Jfunction.

Asan alternative to theexplicit formulation of dependency throughdistributionparameters, itis alsopossibletoconstruct dependenciesimplicitlythrougharithmeticoperators.Forexam- ple,itispossibletorecreatetheexampleaboveusingadditionof stochasticvariablesinsteadoflettingadistributionparameterbe stochastic.Moreprecisely,wehaveauniformvariableon[−1,1]

andanormallydistributedvariablewithlocationatx=0.Adding

theuniformvariabletothenormalvariablecreatesanewnormal variablewithstochasticlocation:

Asbefore, the software automatically sortsthe dependency ordering from the context. Here, since the uniform variable is present as first argument, the software recognises the second argumentasa normal distribution,conditionedontheuniform distribution,andnottheotherwayaround.

Anotherfavorable featurein Chaospyis that multipletrans- formations can be stacked on top of each other. For example, considertheexampleofamultivariatelog-normalrandomvari- able Q withthree dependentcomponents.(Letus ignorefor a momentthefactthatChaospyalreadyofferssuchadistribution.) Tryingtodecomposethisdistributionisaverycumbersometask ifperformedmanually.However,thisprocesscanbedrastically simplified throughvariable transformations,for which Chaospy hasstrongsupport.Alog-normaldistribution,forexample,canbe expressedas

Q=e^ZL+b,

where Zarestandardnormalvariables,andLand bareprede- finedmatrixandvector,respectively.Toimplementthisparticular transformation,weonlyhavetowrite

Theresultingdistributionisfullyfunctionalmultivariatelog- normal,assumingLand bareproperlydefined.

Oneobviousprerequisiteforusingunivariatedistributionsto createconditionalsandmultivariatedistributions,istheavailabil- ityofunivariatedistributions.Sincetheunivariatedistributionis thefundamentalbuildingblock,Chaospyoffersalargecollection of64univariatedistributions.TheyarealllistedinTable1.Thelist alsoshowsthatDakota’ssupportislimitedto11distributions,and Turnshasacollectionof26distributions.

TheChaospysoftwaresupportsinadditioncustomdistributions throughthefunctioncp.constructor.Toillustrateitsuse,con- sider thesimple exampleof a uniformrandomvariableonthe

Table2

Thelistofsupportedcopulasinthevarioussoftwarepackages.

Supportedcopulas Dakota Turns Chaospy

Ali–Mikhail–Haq No Yes Yes

Clayton No Yes Yes

Farlie–Gumbel–Morgenstein No Yes No

Frank No Yes Yes

Gumbel No Yes Yes

Joe No No Yes

Minimum No Yes No

Normal/Nataf Yes Yes Yes

interval[lo,up].Theminimalinputtocreatesuchadistribution is

Here,thetwoprovidedargumentsareacumulativedistribution function(cdf),andaboundaryintervalfunction(bnd),respectively.

Thecp.constructor functionalsotakesseveraloptionalargu- mentstoprovideextrafunctionality.Forexample,theinverseof thecumulativedistributionfunction–thepointpercentilefunction –canbeprovidedthroughtheppfkeyword.Ifthisfunctionisnot provided,thesoftwarewillautomaticallyapproximateitusingthe methoddescribedinSection3.2.

3.4. Copulas

DakotaandTurnsdonotsupporttheRosenblatttransformation appliedtomultivariatedistributionswithdependencies.Instead, thetwopackagesmodeldependenciesusingcopulas[16].Acopula consistsofstochasticallyindependentmultivariatedistributions madedependentusingaparameterizedfunctiong.SincetheRosen- blatttransformationisgeneralpurpose,itispossibletoconstruct anycopuladirectly.However,thiscanquicklybecomeaverycum- bersometasksinceeachcopulamustbedecomposedindividually foreachcombinationofindependentdistributionsandparameter- izationofg.Tosimplifytheuser’sefforts,Chaospyhasdedicated constructorsthatcanreformulateacopulacouplingintoaRosen- blatttransformation.ThisisdonefollowingtheworkofLee[17]and approximatedusingfinitedifferences.Theimplementationisbased ofthesoftwaretoolboxRoseDist[18].Inpractice,thisapproach allowcopulastobedefinedinaRosenblatttransformationsetting.

Forexample,toconstructabivariatenormaldistributionwitha ClaytoncopulainChaospy,wedothefollowing:

A list of supported copulas are listed in Table 2. It shows thatTurnssupports7methods,Chaospy6,whileDakotaoffers1 method.

3.5. Variancereductiontechniques

Asnoted in the beginning of Section 3,by generating sam- ples{Qk}_k∈IK andevaluatingtheresponsefunctionf,itispossible todrawinferenceuponYwithoutknowledgeaboutpY,through MonteCarlosimulation.Unfortunately,thenumberofsamplesK toachieve reasonableaccuracycanoftenbevery high,soiffis assumed tobecomputationallyexpensive, thenumber of sam- plesneededfrequentlymakeMonteCarlosimulationinfeasiblefor practicalapplications.Asawaytomitigatethisproblem,itispos- sibletomodify{Qk}_k∈IK fromtraditionalpseudo-randomsamples,

Table3

Thedifferentsamplingschemesavailable.

Dakota Turns Chaospy

Quasi-MonteCarloscheme

Fauresequence[20] No Yes No

Haltonsequence[21] Yes Yes Yes

Hammersleysequence[22] Yes Yes Yes

Haselgrovesequence[23] No Yes No

Korobovlatice[24] No No Yes

Niederreitersequence[25] No Yes No

Sobolsequence[26] No Yes Yes

Othermethods

Antitheticvariables[1] No No Yes

Importancesampling[1] Yes Yes Yes

LatinHypercubesampling[27] Yes Limited Yes

sothattheaccuracyincreases.Schemesthatselectnon-traditional samplesfor{Qk}_k∈IK toincrease accuracyareknownasvariance reduction techniques.A listof suchtechniques arepresented in Table3,anditshowsthatDakota,TurnsandChaospysupport4, 7,and7variancereductiontechniques,respectively.

Oneofthemorepopularvariance reductiontechniqueisthe quasi-Monte Carlo scheme [1]. The methodconsistsof selecting thesamples {Qk}_k∈IK tobea low-discrepancy sequence instead ofpseudo-randomsamples.Theideaisthatsamplesplacedwith agivendistancefromeachotherincreasethecoverageoverthe samplespace,requiringfewersamplestoreachagivenaccuracy.

Forexample,ifstandardMonteCarlorequires10⁶samplesfora givenaccuracy,quasi-MonteCarlocanoftengetawaywithonly 10³.Notethatthiswouldbreaksomeofthestatisticalproperties ofthesamples[19].

Mostofthetheoryonquasi-MonteCarlomethodsfocuseson generatingsamplesontheunithypercube[0,1]^N.Theoptionto generatesamplesdirectlyontootherdistributionsexists,butis often very limited. To theauthors’ knowledge, the only viable methodforincludingmostquasi-MonteCarlomethodsintothe vastmajorityofnon-standardprobabilitydistributions,isthrough theRosenblatttransformation.SinceChaospyisbuiltaroundthe Rosenblatttransformation,ithasthenovelfeatureofsupporting quasi-MonteCarlomethodsforallprobabilitydistributions.Turns andDakotaonlysupportRosenblatttransformationsforindepen- dentvariablesandtheNormalcopula.

Sometimesthequasi-MonteCarlomethodisinfeasiblebecause theforwardmodelistoocomputationallycostly.Thenextsection describespolynomialchaosexpansions,which oftenrequirefar fewersamplesthanthequasi-MonteCarlomethodforthesame amountofaccuracy.

4. Polynomialchaosexpansions

Polynomialchaosexpansionsrepresentacollectionofmeth- odsthatcanbeconsideredasubsetofpolynomialapproximation methods,butparticularlydesignedforuncertaintyquantification.

Ageneralpolynomialapproximationcanbedefinedas

fˆ(x,t,Q)=

n∈I_N

cn(x,t)˚n(Q) IN={0,...,N}, (11)

where{cn}n∈I_Narecoefficients(oftenknownasFouriercoefficients) and{˚n}_n∈IN arepolynomials.If ˆf isagoodapproximationoff,it ispossibletoeitherinferstatisticalpropertiesof ˆfanalyticallyor throughcheapnumericalcomputationswhere ˆfisusedasasurro- gateforf.

Table4

Methodsforgeneratingexpansionsoforthogonalpolynomials.

OrthogonalizationMethod Dakota Turns Chaospy

Askey–Wilsonscheme[29] Yes Yes Yes

Bertranrecursion[31] No No Yes

Choleskydecomposition[14] No No Yes

DiscretizedStieltjes[32] Yes No Yes

ModifiedChebyshev[32] Yes Yes No

ModifiedGram–Schmidt[32] Yes Yes Yes

A polynomial chaos expansion is defined as a polynomial approximation, as in (11), where the polynomials {˚n}n∈I_N are orthogonalonacustomweightedfunctionspaceLQ:

˚n,˚m =E˚n(Q)˚m(Q)

=

...

˚n(q)˚m(q)pQ(q)dq=0 n=/m. (12)

As a side note, it is worth noting that in parallel with polynomial chaos expansions, there also exists an alternative collocationmethodbasedonmultivariateLagrangepolynomials [28].ThismethodissupportedbyDakotaandChaospy, butnot Turns.

Togenerateapolynomialchaosexpansion,wemustfirstcalcu- latethepolynomials{˚n}_n∈INsuchthattheorthogonalityproperty in(12)issatisfied.ThiswillbethetopicofSection4.1.InSection4.2 weshowhowtoestimatethecoefficients{cn}_n∈IN.Last,inSection 4.7,toolsusedtoquantifyuncertaintyinpolynomialchaosexpan- sionswillbediscussed.

4.1. Orthogonalpolynomialsconstruction

From(12)itfollowsthattheorthogonalitypropertyisnotin generaltransferablebetweendistributions,sinceanewsetofpoly- nomialshastobeconstructedforeachp_Q.Theeasiestapproach toconstructorthogonalpolynomialsistoidentifytheprobability densityp_Q intheso-calledAskey-Wilsonscheme[29].Thepoly- nomialscanthenbepickedfromalist,orbebuiltfromanalytical components.Thecontinuousdistributionssupportedinthescheme include thestandard normal, gamma,beta, and uniformdistri- butionsrespectively through theHermite,Laguerre, Jacobi, and Legendrepolynomialexpansion.Allthethreementionedsoftware toolboxessupporttheseexpansions.

Moving beyondthestandard collectionof theAskey-Wilson scheme,it ispossibletocreatecustom orthogonalpolynomials, both analytically and numerically. Unfortunately, most meth- odsinvolvingfiniteprecisionarithmeticsareill-posed,makinga numericalapproachquiteachallenge[30].Thissectionexplores thevarious approachesfor constructingpolynomialexpansions.

Afulllistof methodsisfoundinTable4.Itshowsthat Dakota, TurnsandChaospysupport4,3and5orthogonalisationmethods, respectively.

Looking beyond an analytical approach, the most popular methodforconstructingorthogonalpolynomialsisthediscretized Stieltjesprocedure[33].Asfarastheauthorsknow,itistheonly trulynumericallystablemethodfororthogonalpolynomialcon- struction.Itisbaseduponone-dimensionalrecursioncoefficients thatareestimatedusingnumericalintegration.Unfortunately,the methodisonlyapplicableinthemultivariatecaseifthecompo- nentsofp_Qarestochasticallyindependent.

Generalizedpolynomialchaosexpansions.Oneapproachtomodel densitieswithstochasticallydependentcomponentsnumerically, istoreformulatetheuncertaintyproblemasasetofindependent

componentsthroughgeneralisedpolynomialchaosexpansion[34].

AsdescribedindetailinSection3.1,aRosenblatttransformation allowsforthemappingbetweenanydomainandtheunithyper- cube[0,1]^D.Withadoubletransformationwecanreformulatethe responsefunctionfas

f(x,t,Q)=f(x,t,T_Q⁻¹(TR(R)))≈fˆ(x,t,R)=

n∈I_N

cn(x,t)˚n(R),

where Risanyrandomvariabledrawnfromp_R,whichforsim- plicity is chosen to consistsof independent components. Also, {˚n}n∈I_N isconstructedtobeorthogonalwithrespecttoL_R,not L_Q. In any case, R is either selected from the Askey-Wilson scheme,or calculated using the discretizedStieltjes procedure.

We remarkthat the accuracy ofthe approximationdeteriorate if the transformation composition T_Q⁻¹◦T_R is not smooth [34].

Dakota, Turns, and Chaospy allsupportgeneralized polynomial chaosexpansionsforindependentstochasticvariablesandtheNor- mal/NatafcopulalistedinTable2.SinceChaospyhastheRosenblatt transformationunderlyingthecomputationalframework,gener- alizedpolynomial chaosexpansions arein fact availablefor all densities.

Thedirectmultivariateapproach.Giventhatboththedensityp_Q hasstochasticallydependentcomponents,andthetransformation compositionT_Q⁻¹◦T_Risnotsmooth,itisstillpossibletogenerate orthogonalpolynomialsnumerically.Asnotedabove,mostmeth- odsarenumericallyunstable,andtheaccuracyintheorthogonality candeterioratewithpolynomialorder,butthemethodscanstillbe useful[14].InTable4,onlyChaospy’simplementationofBertran’s recursion method[31], Choleskydecomposition [35] and mod- ifiedGram-Schmidt orthogonalization[32]supportconstruction of orthogonalpolynomialsfor multivariatedependentdensities directly.

Custompolynomialexpansions.Inthemostextremecases,an automatednumericalmethodisinsufficient.Instead,apolynomial expansionhastobeconstructedmanually.User-definedexpan- sionscan becreatedconveniently, asdemonstratedin thenext exampleinvolvingasecond-orderHermitepolynomialexpansion, orthogonalwithrespecttothenormaldensity[29]:

n∈I6=

1,Q0,Q1,Q₀²−1,Q0Q1,Q₁²−1

TherelevantChaospycodeforcreatingthispolynomialexpansion lookslike

Chaospycontainsacollectionoftoolstomanipulateandcreate polynomials,seeTable5.

Onethingworthnotingisthatpolynomialchaosexpansionssuf- fersfromthecurseofdimensionality:Thenumberoftermsgrows exponentiallywiththenumber ofdimensions [36].Asa result, Chaospydoesnotsupportneitherhighdimensionalnorinfinite dimensionalproblems(randomfields).Oneapproachtoaddress suchproblemswithpolynomialchaosexpansionistofirstreduce thenumberofdimensionthroughtechniqueslikeKarhunen–Loeve expansions[37].Ifsoftwareimplementationsofsuchmethodscan beprovided,theusercaneasilyextendChaospytohighandinfinite dimensionalproblems.

ChaospyincludesoperatorssuchastheexpectationoperatorE. Thisisahelpfultooltoensurethattheconstructedpolynomials areorthogonal,asdefinedin(12).Toverifythattwoelementsin

Table5

Listoftoolsforcreatingandmanipulatingpolynomials.

Function Description

all Testallcoefficientsfornon-zero

any Testanycoefficientsfornon-zero

around Roundtoagivendecimal

asfloat Setcoefficientstypeasfloat

asint Setcoefficienttypeasint

basis Createmonomialbasis

cumprod Cumulativeproduct

cumsum Cumulativesum

cutoff Truncatepolynomialorder

decompose Convertfromseriestosequence

diag Constructorextractdiagonal

differential Differentialoperator

dot Dot-product

flatten Flattenanarray

gradient Gradient(orJacobian)operator

hessian Hessianoperator

inner Innerproduct

mean Average

order Extractpolynomialorder

outer Outerproduct

prod Product

repeat Repeatpolynomials

reshape Reshapeaxes

roll Rollpolynomials

rollaxis Rollaxis

rolldim Rollthedimension

std Empiricalstandarddeviation

substitute Variablesubstitution

sum Sumalonganaxis

swapaxes Interchangetwoaxes

swapdim Swapthedimensions

trace Sumalongthediagonal

transpose Transposethecoefficients

tril Extractlowertriangleofcoefficients

tricu Extractcross-diagonaluppertriangle

var Empiricalvariance

variable Simplepolynomialconstructor

phiare indeedorthogonalunderthestandardbivariatenormal distribution,onewrites

>>> dist = cp.J(cp.Normal(0,1), cp.Normal(0,1))

>>> print cp.E(phi[3]*phi[5], dist) 0.0

Moredetailsofoperatorsusedtoperformuncertaintyanalysis aregiveninSection4.7.

4.2. Calculatingcoefficients

Thereareseveralmethodologiesforestimatingthecoefficients {cn}_n∈IN,typicallycategorizedeitherasnon-intrusiveorintrusive, where non-intrusive means that the computationalprocedures only requires evaluation of f(i.e., software for fcan be reused asa blackbox).Intrusivemethodsneedtoincorporateinforma- tionabouttheunderlyingforwardmodelinthecomputationof thecoefficients. Incaseofforwardmodelsbasedondifferential equations,oneperformsaGalerkinformulationforthecoefficients in probabilityspace,leading effectively toa D-dimensional dif- ferential equation problem in this space [38]. Back et al. [39]

demonstrated that thecomputational costof suchan intrusive Galerkinmethodinsomecaseswashigherthansomenon-intrusive methods.Noneofthethreetoolboxesdiscussedinthispaperhave supportforintrusivemethods.

Withintherealmofnon-intrusivemethods,thereareinprinci- pletwoviablemethodologiesavailable:pseudo-spectralprojection [40] andthepoint collocationmethod[41].Theformer applies anumericalintegrationschemetoestimateFouriercoefficients, whilethelattersolves alinear systemarisingfroma statistical

Table6

Various numerical integration strategies implemented in the three software toolboxes.

Nodeandweightgenerators Dakota Turns Chaospy

Clenshaw-Curtisquadrature[42] Yes No Yes

Cubaturerules[43] Yes No No

Gauss-Legendrequadrature[44] Yes No Yes

Gauss-Pattersonquadrature[45] Yes No Yes

Genz-Keisterquadrature[46] Yes No Yes

Lejaquadrature[47] No No Yes

MonteCarlointegration[1] Yes No Yes

OptimalGaussianquadrature[44] Yes No Yes

regressionformulation.DakotaandChaospysupportbothmethod- ologies,whileTurnsonlysupportspointcollocation.Weshallnow discussthepractical,genericimplementationofthesetwomethods inChaospy.

4.3. Integrationmethods

Thepseudo-spectralprojectionmethodisbasedonastandard leastsquaresminimizationintheweightedfunctionspaceL_Q.Since thepolynomialsareorthogonalinthisspace,theassociatedlin- earsystemisdiagonal,whichallowsaclosed-formexpressionfor theFouriercoefficients.Theexpressioninvolveshigh-dimensional integralsinL_Q.Numericalintegrationisthenrequired,

cn= EY˚n

E˚²n

= 1

E˚²n

...

pQ(q)f(x,t,q)˚n(q)dq

≈ 1

E˚²_n

k∈I_K

w_kp_Q(qk)f(x,t,qk)˚n(qk) IK={0,...,K−1}, (13)

wherewkareweightsand qknodesinaquadraturescheme.Note thatfisonlyevaluatedforthenodes qk,andtheseevaluationscan bemadeonce.Thereafter,onecanexperimentwiththepolynomial ordersinceanyc_ndependsonthesameevaluationsoff.

Table 6 shows the various quadrature schemes offered by Dakota and Chaospy (recall that Turns does not sup- port pseudo-spectral projection). All techniques for generating nodes and weights in Chaospy are available through the cp.generatequadraturefunction.Supposewewanttogenerate optimalGaussianquadraturenodesforthenormaldistribution.We thenwrite

Mostquadratureschemesaredesignedforunivariateproblems.

Toextendaunivariateschemetothemultivariatecase,integra- tionrulesalongeachaxiscanbecombinedusingatensorproduct.

Unfortunately,suchaproductsuffersfromthecurseofdimension- ality andbecomes a verycostlyintegrationprocedurefor large D.In higher-dimensionalproblemsonecanreplacethefullten- sorproductbyaSmolyaksparsegrid[48].Themethodworksby takingmultiplelowerordertensorproductrulesandjoiningthem together.Iftheruleisnested,i.e.,thesamesamplesfoundatalow orderarealsoincludedathigherorder,thenumberofevaluations canbefurtherreduced.Anotherfeatureistoaddanisotropysuch thatsomedimensionsaresampledmorethanothers[49].Inaddi- tiontothetensorproductrules,thereareafewnativemultivariate cubaturerulesthat allowforlow ordermultivariateintegration [43].BothDakotaandChaospyalsosupporttheSmolyaksparse gridandanisotropy.

Chaospy hassupport for construction of custom integration rulesdefinedbytheuser.Thecp.rulegeneratorfunctioncanbe usedtojoinalistofunivariaterulesusingtensorgridorSmolyak sparsegrid.Forexample,considerthetrapezoidrule:

Thecp.rule generatorfunctiontakespositionalarguments, eachrepresentingaunivariaterule.Togeneratearuleforthemul- tivariatecase,withthesameone-dimensionalrulealongtwoaxes, wedothefollowing:

Softwareforconstructingandexecutingageneral-purposeinte- grationschemeisusefulforseveralcomputationalcomponentsin uncertaintyquantification.Forexample,inSection4.1whencon- structingorthogonalpolynomialsusingrawstatisticalmoments,or calculatingdiscretizedStieltjes’recurrencecoefficients,numerical integrationisrelevant.LiketheppffunctionnotedinSection3.3, themomentsandrecurrencecoefficientscanbeaddeddirectlyinto eachdistribution.However,whenthesearenotavailable,Chaospy will automaticallyestimate missing information by quadrature rules, using the cp.generatequadrature function described above.

TocomputetheFouriercoefficientsandthepolynomialchaos expansion,weusethecp.fitquadraturefunction.Ittakesfour arguments:thesetoforthogonalpolynomials,quadraturenodes, quadratureweights,andtheuser’sfunctionforevaluatingthefor- wardmodel(tobeexecutedatthequadraturenodes).Notethatin thecaseofthediscretizedStieltjesmethoddiscussedinSection4.1, thenominatorE˚²_nin(13)canbecalculatedmoreaccuratelyusing recurrencecoefficients[32].Specialnumericalfeatureslikethiscan beaddedbyincludingoptionalargumentsincp.fitquadrature.

4.4. Pointcollocation

Theothernon-intrusiveapproachtoestimatethecoefficients {ck}_k∈IKisthepointcollocationmethod.Onewayofformulatingthe methodistorequirethepolynomialexpansiontoequalthemodel evaluationsatasetofcollocationnodes{Qk}_k∈IK,resultinginan over-determinedsetoflinearequationsfortheFouriercoefficients:

⎡

⎢ ⎢

⎣

0(q0) ··· N(q0) ..

. ...

0(q_K−1) ··· N(q_K−1)

⎤

⎥ ⎥

⎦

⎡

⎢ ⎢

⎣

c0

.. . cN

⎤

⎥ ⎥

⎦

⎡

⎢ ⎢

⎣

f(q0) .. . f(q_K−1)

⎤

⎥ ⎥

⎦

Unlike pseudo spectral projection, thelocations ofthe colloca- tionnodesarenotrequiredtofollowanyintegrationrule.Hosder [41]showedthat thesolutionusingHammersleysamplesfrom quasi-MonteCarlosamplesresultedin morestableresultsthan usingconventionalpseudo-randomsamples.Inotherwords,well placedcollocationnodesmightincreasetheaccuracy.InChaospy thesecollocationnodescanbeselectedfromintegrationrulesor

Table7

ListofsupportedregressionmethodsforestimatingFouriercoefficients.

Regressionschemes Dakota Turns Chaospy

Basispursuit[52] Yes No No

Bayesianauto.relevancedetermination[53] No No Yes

Bayesianridge[54] No No Yes

Elasticnet[55] Yes No Yes

Forwardstagewise[56] No Yes No

Leastabsoluteshrinkage&Selection[51] Yes Yes Yes Leastangle&ShrinkagewithAIC/BIC[57] No No Yes

Leastsquaresminimization Yes Yes Yes

Orthogonalmatchingpursuit[58] Yes No Yes

Singularvaluedecomposition No Yes No

Tikhonovregularization[50] No No Yes

frompseudo-randomsamplesfromMonteCarlosimulation,asdis- cussedinSection3.5.Inaddition,thesoftwareacceptsuserdefined strategiesforchoosingthesamplingpoints.Turnsalsoallowsfor user-definedpoints,whileDakotahasitspredefinedstrategies.

Theobviouswaytosolvetheover-determinedsystemin(14)is touseleastsquaresminimization,whichresemblesthestandard statisticallinearregressionapproachoffittingapolynomialtoaset ofdatapoints.However,fromanumericalpointofview,thismight notbethebeststrategy.Ifthenumericalstabilityofthesolutionis low,itmightbeprudenttouseTikhonovregularization[50],orif theproblemissolargethatthenumberofcoefficientsisveryhigh,it mightbeusefultoforcesomeofthecoefficientstobezerothrough leastangleregression[51].Beingabletorunandcomparealter- nativemethodsisimportantinmanyproblemstoseeifnumerical stabilityisapotentialproblem.Table7liststheregressionmethods offeredbyDakota,Turns,andChaospy.

Generatingapolynomialchaosexpansionusinglinearregres- sionisdoneusingChaospy’scp.fitregressionfunction.Ittakes thesameargumentsascp.fitquadrature,exceptthatquadra- tureweightsareomitted,andoptionalargumentsdefinetherule usedtooptimize(14).

4.5. Modelevaluations

Irrespectivelyofthemethodusedtoestimatethecoefficientsc_k, theuserisleftwiththejobtoevaluatetheforwardmodel(response function)f,whichisnormallybyfarthemostcomputing-intensive partinuncertaintyquantification.Chaospydoesnotimposeany restrictiononthesimulationcodeusedtocomputetheforward model.Theonlyrequirementisthattheusercanprovideanarray ofvaluesoffatthequadratureorcollocationnodes.Chaospyusers willusuallywrapanycomplexsimulationcodeforfinaPython functionf(q),whereqisanodeinprobabilityspace(i.e.,qcontains valuesoftheuncertainparametersintheproblem).Forexample, forpseudo-spectralprojection,samplesoffcanbecreatedas

orperhapsdoneinparalleliffistimeconsumingtoevaluate:

Theevaluationofallthefvaluescanalsobedoneinparallelwith MPIinadistributedwayonaclusterusingthePythonmodulelike mpi4py.BothDakotaandTurnssupportparallelevaluationoffval- ues,butthefeatureisembededintothecode,potentiallylimiting thecustomizationoptionsoftheparallelization.

4.6. Extensionofpolynomialexpansions

Thereismuchliteraturethatextendsonthetheoryofpolyno- mialchaosexpansion[36].Forexample,Isukapallishowedthatthe

accuracyofapolynomialexpansioncouldbeincreasedbyusing partialderivativesofthemodelresponse[59].Thistheoryisonly directlysupportedbyDakota.InTurnsandChaospythesupportis indirectbyallowingtheusertoaddthefeaturemanually.

Tobeabletoincorporatepartialderivativesoftheresponse,the partialderivativeofthepolynomialexpansionmustbeavailableas well.InbothTurnsandChaospy,thederivativeofapolynomialcan begeneratedeasily.Thisderivativecanthenbeaddedtotheexpan- sion,allowingustoincorporateIsukapalli’stheoryinpractice.This isjustanexampleonhowmanipulationofthepolynomialexpan- sionsandmodelapproximationscanovercomethelackofsupport foraparticularfeaturefromtheliterature.

To be able to support many current and possible future extensions of polynomial chaos, a large collection of tools for manipulatingpolynomialsmustbeavailable.InDakota,nosuch toolsexistfromauserperspective.InTurns,thereissupportfor somearithmeticoperatorsinadditiontothederivative.InChaospy, however,thepolynomialgeneratedforthemodelresponseisofthe sametypeasthepolynomialsgeneratedinSections4.1and4.2,and therichsetofmanipulationsofpolynomialsisthenavailablefor ˆf aswell.

Beyondtheanalyticaltoolsforstatisticalanalysisof ˆf,either fromthetoolboxorcustomonesbytheuser,therearemanystatis- ticalmetricsthatcannoteasilybeexpressedassimpleclosed-form formulas. Such metrics include confidence intervals, sensitivity indices,p-valuesinhypothesistesting,tomentionafew.Inthose scenarios,itmakessensetoperformasecondaryuncertaintyanaly- sisthroughMonteCarlosimulation.Evaluatingtheapproximation ˆf isnormallycomputationallymuchcheaperthanevaluatingthe fullforwardmodelf,thusallowingalargenumberofMonteCarlo sampleswithinacheapcomputationalbudget.Thistypeofsec- ondarysimulationsaredoneautomaticallyinthebackgroundin DakotaandTurns,whileChaospydoesnotfeatureautomatedtools forsecondaryMonteCarlosimulation.Instead,Chaospyallowsfor simpleandcomputationallycheapgenerationofpseudo-random samples,asdescribedinSection3.5,suchthattheusercaneasily puttogetheratailoredMonteCarlosimulationtomeettheneeds athand.WithinafewlinesofPythoncode,thesamplescanbe analyzedwiththestandardNumpyandtheScipylibraries[60]or withmorespecializedstatisticallibrarieslikePandas[9],Scikit- learn[61],Scikit-statsmodel[62],andPython’sinterfacetotherich Renvironmentforstatisticalcomputing.Forexample,forthespe- cific ˆffunctionillustratedabove,thefollowingcodecomputesa90 percentconfidenceinterval,basedon10⁵pseudo-randomsamples andNumpy’sfunctionalityforfindingpercentilesindiscretedata:

Sincethetypeofstatisticalanalysisof ˆfoftenstronglydepends onthephysicalproblemathand,webelieve thattheabilityto quicklycomposecustomsolutionsbyputtingtogetherbasicbuild- ingblocksisveryusefulinuncertaintyquantification.Thisisyet anotherexampleoftheneedforapackagewithastrongfocuson easycustomization.

4.7. Descriptivetools

Thelast step inuncertainty quantification basedonpolyno- mialchaosexpansionsistoquantifytheuncertainty.Inpolynomial chaosexpansionthisisdonebyusingtheuncertaintyinthemodel approximation f approx as a substiute for the uncertainty in themodelf.Forthemostpopularstatisticalmetrics,likemean,

Table8

Listofcommonstatisticaloperatorsthatcanbeusedforanalyticalevaluationof polynomials.

Method Dakota Turns Chaospy

Covariance/correlation Yes Yes Yes

Expectedvalue Yes Yes Yes

Conditionalexpectation No No Yes

Kurtosis Yes Yes Yes

Sensitivityindex Yes Yes Yes

Skewness Yes Yes Yes

Variance Yes Yes Yes

variance,correlation,apolynomialchaosexpansionallowsforana- lyticalanalysis,whichiseasytocalculateandhashighaccuracy.

Thispropertyisreflectedinallthethreetoolboxes.Tocalculate the expected value,variance and correlation of a simple (here univariate)polynomialapproximationfapprox,withanormally distributed0variable,wecanwithChaospywrite

AlistofsupportedanalyticalmetricsislistedinTable8.

5. Conclusionandfurtherwork

Untilnowtherehaveonlybeenafewrealsoftwarealternatives forimplementingnon-intrusivepolynomialchaosexpansions.Two ofthemorepopularimplementations,DakotaandTurns,areboth high-qualitysoftwarethatcanbeappliedtoalargearrayofprob- lems.Thepresentpaperhasintroducedanewalternative:Chaospy.

Its aimis tobean experimentalfoundryfor scientists. Besides featuringa vastlibraryof state-of-the-arttools,Chaospyallows forahighdegreeofcustomizationinauser-friendlyway.Within afew linesofhigh-levelPythoncode,theusercanplayaround withcustomdistributions,custompolynomials, customintegra- tionschemes, customsamplingschemes,and customstatistical analysisoftheresult.Throughoutthetextwehavecomparedthe built-infunctionalityofthethreepackages,andChaospydovery wellinthiscomparison,whichissummarizedinTable9.Butthe primaryadvantageofthepackageisthestrongemphasisonoffer- ingwell-designedsoftwarebuildingblocks,withahighabstraction level,thatcaneasilybecombinedtocreatetailoreduncertainty quantificationalgorithmsfornewproblems.

Althoughtheprimaryaimofthesoftwareistoconstructpoly- nomialchaosexpansions,thesoftwareis alsoastate-of-the-art toolboxforperformingMonteCarlosimulation,eitherdirectlyon theforwardmodelorincombinationwithpolynomialchaosexpan- sions.Variancereductiontechniquesareincludedtospeedupthe

Table9

AsummaryofthevariousfeaturesinDakota,TurnsandChaospy.

Feature Dakota Turns Chaospy

Distributions 11 26 64

Copulas 1 7 6

Samplingschemes 4 7.5 7

Orthogonalpolynomialschemes 4 3 5

Numericalintegrationstrategies 7 0 7

Regressionmethods 5 4 8

Analyticalmetrics 6 6 7