Hi Pamphil,
You were almost right: the AdaptiveStieltjesAlgorithm is very close to what you
are looking for, but not exactly what you need. It is the algorithmic part of
the factory of orthonormal polynomials, the class you have to use is
StandardDistributionPolynomialFactory, ie a factory (=able to build something)
and not an algorithm (=something able to compute something). You have all the
details here:
http://openturns.github.io/openturns/master/user_manual/_generated/openturns.StandardDistributionPolynomialFactory.html
I agree on the fact that the difference is quite subtle, as it can be seen by
comparing the API of the two classes. The distinction was made at a time were
several algorithms were competing for the task (GramSchmidtAlgorithm,
ChebychevAlgorithm) but in fact the AdaptiveStieltjesAlgorithm proved to be
much more accurate and reliable than the other algorithms, and now it is the
only orthonormalization algorithm available.
Another subtle trick is the following.
If you create a basis this way:
basis = ot.StandardDistributionPolynomialFactory(dist)
you will get the basis associated to the *standard representative* distribution
in the parametric family to which dist belongs. It means the distribution with
zero mean and unit variance, or with support equals to [-1, 1], or dist itself
if no affine transformation is able to reduce the number of parameters of the
distribution.
It is managed automatically within the FunctionalChaosAlgorithm, but can be
disturbing if you do things by hand.
If you create a basis this way:
basis =
ot.StandardDistributionPolynomialFactory(ot.AdaptiveStieltjesAlgorithm(dist))
then the distribution is preserved, and you get the orthonormal polynomials
corresponding to dist. Be aware of the fact that the algorithm may have hard
time to build the polynomials if your distribution is far away from its
standard representative, as it may involves the computation of recurrence
coefficients with a much wider range of variation. The benefit is that the
orthonormality measure is exactly your distribution, assuming that its copula
is the independent one, so you don't have to introduce a marginal
transformation between both measures.
Some additional remarks:
+ it looks like dist has dimension>1, as you extract its marginal distributions
later on. AdaptiveStieltjesAlgorithm and StandardDistributionPolynomialFactory
only work with 1D distributions (it is not checked by the library, my shame).
What you have to do is:
basis =
ot.OrthogonalProductPolynomialFactory([ot.StandardDistributionPolynomialFactory(ot.AdaptiveStieltjesAlgorithm(dist.getMarginal(i)))
for i in range(dist.getDimension())])
Quite a long line, I know...
It will build a multivariate polynomial basis orthonormal with respect to the
product distribution (ie with independent copula) sharing the same 1D marginal
distributions as dist.
After that, everything will work as expected and you will NOT have to build the
transformation (if you build it it will coincide with the identity function).
If you encounter performance issues (the polynomials of high degrees take ages
to be built as in http://trac.openturns.org/ticket/885, or there is an
overflow, or the numerical precision is bad) then use:
basis =
ot.OrthogonalProductPolynomialFactory([ot.StandardDistributionPolynomialFactory(dist.getMarginal(i))
for i in range(dist.getDimension())])
and build the transformation the way you do it.
+ if you use the FunctionalChaosAlgorithm class by providing an input sample
and an output sample, you also have to provide the weights of the input sample
EVEN IF the experiment given in the projection strategy would allow to
recompute them. It is because the fact that you provide the input sample
overwrite the weighted experiment of the projection stratey by a
FixedExperiment doe.
I attached two complete examples: one using the exact marginal distributions
and the other using the standard representatives.
Best regards
Régis
________________________________
De : roy <[email protected]>
À : regis lebrun <[email protected]>
Cc : users <[email protected]>
Envoyé le : Vendredi 6 octobre 2017 14h22
Objet : Re: [ot-users] Sample transformation
Hi Regis,
Thank you for this detailed answer.
- I am using the latest release from conda (OT 1.9, python 3.6.2, latest numpy,
etc.) ,
- For the sample, I need it to generate externally the output (cost code that
cannot be integrated into OT as model),
- I have to convert ot.Sample into np.array because it is then used by other
functions to create the simulations, etc.
If I understood correctly, I can create the projection strategy using this
snippet:
basis = ot.AdaptiveStieltjesAlgorithm(dist)
measure = basis.getMeasure()
quad = ot.Indices(in_dim)
for i in range(in_dim):
quad[i] = degree + 1
comp_dist = ot.GaussProductExperiment(measure, quad)
proj_strategy = ot.IntegrationStrategy(comp_dist)
inv_trans =
ot.Function(ot.MarginalTransformationEvaluation([measure.getMarginal(i) for i
in range(in_dim)], distributions))
sample = np.array(inv_trans(comp_dist.generate()))
It seems to work. Except that the basis does not work with
ot.FixedStrategy(basis, dim_basis). I get a non implemented method error.
After I get the sample and the corresponding output, what is the way to go?
Which arguments do I need to use for the
ot.FunctionalChaosAlgorithm?
I am comparing the Q2 and on Ishigami and I was only able to get correct
results using:
pc_algo = ot.FunctionalChaosAlgorithm(sample, output, dist, trunc_strategy)
But for least square strategy I had to use this:
pc_algo = ot.FunctionalChaosAlgorithm(sample, output)
Is it normal?
Pamphile ROY
Chercheur doctorant en Quantification d’Incertitudes
CERFACS - Toulouse (31) - France
+33 (0) 5 61 19 31 57
+33 (0) 7 86 43 24 22
Le 5 oct. 2017 à 15:40, regis lebrun <[email protected]> a
écrit :
>
>Hi Pamphile,
>
>
>1) The problem:
>The problem you get is due to the fact that in your version of OpenTURNS (1.7
>I suppose), the GaussProductExperiment class has a different way to handle the
>input distribution than the other WeightedExperiment classes: it generates the
>quadrature rule of the *standard representatives* of the marginal
>distributions instead of the marginal distributions. It does not change the
>rate of convergence of the PCE algorithm and allows to use specific algorithms
>for distributions with known orthonormal polynomials. It is not explained in
>the documentation and if you ask the doe for its distribution it will give you
>the initial distribution instead of the standardized one.
>
>2) The mathematical background:
>The generation of quadrature rules for arbitrary 1D distributions is a badly
>conditioned problem. Even if the quadrature rule is well-defined (existence of
>moments of any order, distribution characterized by these moments), the
>application that maps the recurrence coefficients of the orthogonal
>polynomials to their value can have a very large condition number. As a
>result, the adaptive integration used to compute the recurrence coefficients
>of order n, based on the values of the polynomials of degree n-1 and n-2, can
>lead to wrong values and all the process falls down.
>
>3) The current state of the software:
>Since version 1.8, OpenTURNS no more generates the quadrature rule of the
>standard representatives, but the quadrature rule of the actual marginal
>distributions. The AdaptiveStieltjesAlgorithm class, introduced in release
>1.8, is much more robust than the previous orthonormalization algorithms and
>is able to handle even stiff problems. There are still difficult situations
>(distributions with discontinuous PDF inside of the range, fixed in OT 1.9, or
>really badly conditioned distributions, hopefully fixed when ticket#885 will
>be solved) but most usual situations are under control even with marginal
>degrees of order 20.
>
>4) The (probable) bug in your code and the way to solve it
>You must be aware of the fact that the distribution you put into your
>WeightedExperiment object will be superseded by the distribution corresponding
>to your OrthogonalBasisFactory inside of the FunctionalChaosAlgorithm. If you
>need to have the input sample before to run the functional chaos algorithm,
>then you have to build your transformation by hand. Assuming that you already
>defined your projection basis called 'myBasis', your marginal integration
>degrees 'myDegrees' and your marginal distributions 'myMarginals', you have to
>write (in OT 1.7):
>
># Here the explicit cast into a NumericalMathFunction is to be able to
>evaluate the transformation over a sample
>myTransformation =
>ot.NumericalMathFunction(ot.MarginalTransformationEvaluation([myBasis.getDistribution().getMarginal(i)
> for i in range(dimension), myMarginals))
>sample = myTransformation(ot.GaussProductExperiment(myBasis.getDistribution(),
>myDegrees).generate())
>
>
>You should avoid to cast OT objects into np objects as much as possible, and
>if you cannot avoid these casts you should do them only in the sections where
>they are needed. They can be expansive for large objects, and if the sample
>you get from generate() is used only as an argument of a
>NumericalMathFunction, then it will be converted back into a NumericalSample!
>
>Best regards
>
>Régis
>________________________________
>De : roy <[email protected]>
>À : users <[email protected]>
>Envoyé le : Jeudi 5 octobre 2017 11h13
>Objet : [ot-users] Sample transformation
>
>
>
>Hi,
>
>I am facing consistency concerns in the API regarding distributions and
>sampling.
>
>The initial goal was to get the sampling for Polynomial Chaos as I must not
>use the model variable.
>So for least square strategy I do something like this:
>
>proj_strategy = ot.LeastSquaresStrategy(montecarlo_design)
>sample = np.array(proj_strategy.getExperiment().generate())
>
>sample is correct as the bounds of each feature lie in the corresponding
>ranges.
>
>But now if I want to use IntegrationStrategy:
>
>ot.IntegrationStrategy(ot.GaussProductExperiment(dists, list))
>sample = np.array(proj_strategy.getExperiment().generate())
>
>sample’s outputs lie between [-1, 1] which does not corresponds to the
>distribution I have initially.
>
>So I used the conversion class but it does not work well with
>GaussProductExperiment as it requires [0, 1] instead of [-1, 1].
>
>Thus I use this hack:
>
># Convert from [-1, 1] -> input distributions
>marg_inv_transf = ot.MarginalTransformationEvaluation(distributions, 1)
>sample = (proj_strategy.getExperiment().generate() + 1) / 2.
>
>
>Is it normal that the distribution classes are not returning in the same
>intervals?
>
>
>Thanks for your support!
>
>
>Pamphile ROY
>Chercheur doctorant en Quantification d’Incertitudes
>CERFACS - Toulouse (31) - France
>+33 (0) 5 61 19 31 57
>+33 (0) 7 86 43 24 22
>
>
>_______________________________________________
>OpenTURNS users mailing list
>[email protected]
>http://openturns.org/mailman/listinfo/users
>import openturns as ot
# The model
model = ot.SymbolicFunction(["x", "y"], ["sin(x) / (1.0 + y * y)"])
# The distribution and the associated multivariate polynomial basis
marginals = [ot.Exponential(2.0), ot.Normal(-1.0, 0.5)]
dimension = len(marginals)
distribution = ot.ComposedDistribution(marginals)
enumerateFunction = ot.EnumerateFunction(dimension)
basis = ot.OrthogonalProductPolynomialFactory([ot.StandardDistributionPolynomialFactory(ot.AdaptiveStieltjesAlgorithm(margin)) for margin in marginals], enumerateFunction)
measure = basis.getMeasure()
# The Functional Chaos details...
deg = 20
degrees = [deg+1]*dimension
projectionStrategy = ot.IntegrationStrategy(ot.GaussProductExperiment(measure, degrees))
adaptiveStrategy = ot.FixedStrategy(basis, enumerateFunction.getStrataCumulatedCardinal(deg))
# Generate the input sample WITH THE ASSOCIATED WEIGTHS
inputSample, weights = projectionStrategy.getExperiment().generateWithWeights()
# Evaluate the model
outputSample = model(inputSample)
# If you have access to the full model
# algo = ot.FunctionalChaosAlgorithm(model, distribution, adaptiveStrategy, projectionStrategy)
# If you forget to give the weights
# algo = ot.FunctionalChaosAlgorithm(inputSample, outputSample, distribution, adaptiveStrategy, projectionStrategy)
# If you do things the right way
algo = ot.FunctionalChaosAlgorithm(inputSample, weights, outputSample, distribution, adaptiveStrategy, projectionStrategy)
ot.Log.Show(ot.Log.INFO)
algo.run()
result = algo.getResult()
# Some checks
metaModel = result.getMetaModel()
lower = distribution.computeQuantile(0.001)
upper = distribution.computeQuantile(0.999)
graph = model.draw(lower, upper)
graph.setColors(["red"])
graph2 = metaModel.draw(lower, upper)
graph2.setColors(["blue"])
graph.add(graph2)
graph.setLegendPosition("")
ot.Show(graph)
import openturns as ot
# The model
model = ot.SymbolicFunction(["x", "y"], ["sin(x) / (1.0 + y * y)"])
# The distribution and the associated multivariate polynomial basis
marginals = [ot.Exponential(2.0), ot.Normal(-1.0, 0.5)]
dimension = len(marginals)
distribution = ot.ComposedDistribution(marginals)
enumerateFunction = ot.EnumerateFunction(dimension)
basis = ot.OrthogonalProductPolynomialFactory([ot.StandardDistributionPolynomialFactory(margin) for margin in marginals], enumerateFunction)
measure = basis.getMeasure()
transformation = ot.Function(ot.MarginalTransformationEvaluation([measure.getMarginal(i) for i in range(dimension)], marginals, False))
# The Functional Chaos details...
deg = 20
degrees = [deg+1]*dimension
projectionStrategy = ot.IntegrationStrategy(ot.GaussProductExperiment(measure, degrees))
adaptiveStrategy = ot.FixedStrategy(basis, enumerateFunction.getStrataCumulatedCardinal(deg))
# Generate the input sample WITH THE ASSOCIATED WEIGTHS
standardInputSample, weights = projectionStrategy.getExperiment().generateWithWeights()
inputSample = transformation(standardInputSample)
# Evaluate the model
outputSample = model(inputSample)
# If you have access to the full model
# algo = ot.FunctionalChaosAlgorithm(model, distribution, adaptiveStrategy, projectionStrategy)
# If you forget to give the weights
# algo = ot.FunctionalChaosAlgorithm(inputSample, outputSample, distribution, adaptiveStrategy, projectionStrategy)
# If you do things the right way
algo = ot.FunctionalChaosAlgorithm(inputSample, weights, outputSample, distribution, adaptiveStrategy, projectionStrategy)
ot.Log.Show(ot.Log.INFO)
algo.run()
result = algo.getResult()
# Some checks
metaModel = result.getMetaModel()
lower = distribution.computeQuantile(0.001)
upper = distribution.computeQuantile(0.999)
graph = model.draw(lower, upper)
graph.setColors(["red"])
graph2 = metaModel.draw(lower, upper)
graph2.setColors(["blue"])
graph.add(graph2)
graph.setLegendPosition("")
ot.Show(graph)
_______________________________________________
OpenTURNS users mailing list
[email protected]
http://openturns.org/mailman/listinfo/users
