# Estimate moments from Taylor expansions¶

In this example we are going to estimate mean and standard deviation of an output variable of interest thanks to the Taylor variance decomposition method of order 1 or 2.

## Model definition¶

```import openturns as ot
import openturns.viewer as viewer
import numpy as np
from matplotlib import pylab as plt

ot.Log.Show(ot.Log.NONE)
```

Create a composite random vector

```ot.RandomGenerator.SetSeed(0)
input_names = ["x1", "x2", "x3", "x4"]
myFunc = ot.SymbolicFunction(input_names, ["cos(x2*x2+x4)/(x1*x1+1+x3^4)"])
R = ot.CorrelationMatrix(4)
for i in range(4):
R[i, i - 1] = 0.25
distribution = ot.Normal([0.2] * 4, [0.1, 0.2, 0.3, 0.4], R)
distribution.setDescription(input_names)
# We create a distribution-based RandomVector
X = ot.RandomVector(distribution)
# We create a composite RandomVector Y from X and myFunc
Y = ot.CompositeRandomVector(myFunc, X)
```

## Taylor expansion based estimation of the moments¶

We create a Taylor expansion method to approximate moments

```taylor = ot.TaylorExpansionMoments(Y)
```

## Analysis of the results¶

get 1st order mean

```print(taylor.getMeanFirstOrder())
```
```[0.932544]
```

get 2nd order mean

```print(taylor.getMeanSecondOrder())
```
```[0.820295]
```

get covariance

```print(taylor.getCovariance())
```
```[[ 0.0124546 ]]
```

draw importance factors

```taylor.getImportanceFactors()
```
class=PointWithDescription name=Unnamed dimension=4 description=[x1,x2,x3,x4] values=[0.181718,0.0430356,0.0248297,0.750417]

draw importance factors

```graph = taylor.drawImportanceFactors()
view = viewer.View(graph)
```

Get the value of the output at the mean point

```taylor.getValueAtMean()
```
class=Point name=Unnamed dimension=1 values=[0.932544]

Get the gradient value of the output at the mean point

```taylor.getGradientAtMean()
```

[[ -0.35812 ]
[ -0.0912837 ]
[ -0.0286496 ]
[ -0.228209 ]]

Get the hessian value of the output at the mean point

```taylor.getHessianAtMean()
plt.show()
```

When no gradient and/or functions are provided for the model, a finite difference approach is relied on automatically. However, it is possible to manually specify the characteristic of the considered difference steps.

```def myPythonFunction(X):
x1, x2, x3, x4 = X
return [np.cos(x2 * x2 + x4) / (x1 * x1 + 1.0 + x3**4)]

myFunc = ot.PythonFunction(4, 1, myPythonFunction)
```

For instance, a user-defined constant step value can be considered

```gradEpsilon = [1e-8] * 4
hessianEpsilon = [1e-7] * 4
hessianStep = ot.ConstantStep(hessianEpsilon)  # Constant Hessian step
)
myFunc.setHessian(
ot.CenteredFiniteDifferenceHessian(hessianStep, myFunc.getEvaluation())
)
```

Alternatively, we can consider a finite difference step value which depends on the location in the input space by relying on the BlendedStep class:

```gradEpsilon = [1e-8] * 4
hessianEpsilon = [1e-7] * 4
hessianStep = ot.BlendedStep(hessianEpsilon)  # Constant Hessian step
```Y = ot.CompositeRandomVector(myFunc, X)