Metropolis Hastings algorithm

class ffpack.rpm.metropolisHastings.AuModifiedMHSampler(initialVal=None, targetPdf=None, proposalCSampler=None, sampleDomain=None, randomSeed=None, **sdKwargs)

Modified Metropolis-Hastings sampler based on Au and Beck algorithm [Au2001].

References

Au2001: Au, S.K. and Beck, J.L., 2001. Estimation of small failure probabilities in high dimensions by subset simulation. Probabilistic engineering mechanics, 16(4), pp.263-277.

__init__(initialVal=None, targetPdf=None, proposalCSampler=None, sampleDomain=None, randomSeed=None, **sdKwargs)

Initialize the Au modified Metropolis-Hastings sampler

Parameters

initialVal (array_like) – Initial observed data point.
targetPdf (function list) – Target probability density function list. Each element targetPdf[ i ] in the list is a callable function referring the independent marginal. targetPdf[ i ] takes one input parameter and return the corresponding probability. It will be called as targetPdf[ i ]( X[ i ] ) where X is a list in which the element is same type as initialVal[ i ], and a scalar value of probability should be returned by targetPdf[ i ]( X[ i ] ).
proposalCSampler (function list) – Proposal conditional sampler list (i.e., transition kernel list). Each element proposalCSampler[ i ] in the list is a callable function referring a sampler that will return a sample for the given observed data point. A usual choice is to let proposalCSampler[ i ] be a Gaussian/normal distribution centered at the observed data point. It will be called as proposalCSampler[ i ]( X[ i ] ) where X is a list in which each element is the same type as initialVal[ i ], and a sample with the same type of initialVal[ i ] should be returned.
sampleDomain (function, optional) – Sample domain function. sampleDomain is a function to determine if a sample is in the sample domain. For example, if the sample doamin is [ 0, inf ] and the sample is -2, the sample will be rejected. For the sampling on field of real numbers, it should return True regardless of the sample value. It called as sampleDomain( cur, nxt, **sdKwargs ) where cur, nxt are lists in which each element is the same type as initivalVal[ i ], and a boolean value should be returned.
randomSeed (integer, optional) – Random seed. If randomSeed is none or is not an integer, the random seed in global config will be used.

Raises

ValueError – If initialVal, targetPdf, or proposalCSampler is None. If dims of initialVal, targetPdf, and proposalCSampler are not equal. If targetPdf returns negative value.

Examples

>>> from ffpack.rpm import AuModifiedMHSampler
>>> initialValList = [ 1.0, 1.0 ]
>>> targetPdf = [ lambda x : 0 if x < 0 else np.exp( -x ),
...               lambda x : 0 if x < 0 else np.exp( -x ) ]
>>> proposalCSampler = [ lambda x : np.random.normal( x, 1 ),
...                      lambda x : np.random.normal( x, 1 ) ]
>>> auMMHSampler = AuModifiedMHSampler( initialVal, targetPdf,
...                                     proposalCsampler )

getSample()

Get a sample.

Returns: rst – Data point sample.
Return type: array_like

Examples

>>> sample = auMMHSampler.getSample()

class ffpack.rpm.metropolisHastings.MetropolisHastingsSampler(initialVal=None, targetPdf=None, proposalCSampler=None, sampleDomain=None, randomSeed=None, **sdKwargs)

Metropolis-Hastings sampler to sample data [Bourinet2018].

References

Bourinet2018: Bourinet, J.M., 2018. Reliability analysis and optimal design under uncertainty-Focus on adaptive surrogate-based approaches (Doctoral dissertation, Université Clermont Auvergne).

__init__(initialVal=None, targetPdf=None, proposalCSampler=None, sampleDomain=None, randomSeed=None, **sdKwargs)

Initialize the Metropolis-Hastings sampler

Parameters

initialVal (scalar or array_like) – Initial observed data point.
targetPdf (function) – Target probability density function or target distribution function. targetPdf takes one input parameter and return the corresponding probability. It will be called as targetPdf( X ) where X is the same type as initialVal, and a scalar value of probability should be returned.
proposalCSampler (function) – Proposal conditional sampler (i.e., transition kernel). proposalCSampler is a sampler that will return a sample for the given observed data point. A usual choice is to let proposalCSampler be a Gaussian/normal distribution centered at the observed data point. It will be called as proposalCSampler( X ) where X is the same type as initialVal, and a sample with the same type of initialVal should be returned.
sampleDomain (function, optional) – Sample domain function. sampleDomain is a function to determine if a sample is in the sample domain. For example, if the sample doamin is [ 0, inf ] and the sample is -2, the sample will be rejected. For the sampling on field of real numbers, it should return True regardless of the sample value. It called as sampleDomain( cur, nxt, **sdKwargs ) where cur, nxt are the same type as initivalVal, and a boolean value should be returned.
randomSeed (integer, optional) – Random seed. If randomSeed is none or is not an integer, the random seed in global config will be used.

Raises

ValueError – If initialVal, targetPdf, or proposalCSampler is None. If targetPdf returns negative value.

Examples

>>> from ffpack.rpm import MetropolisHastingsSampler
>>> initialVal = 1.0
>>> targetPdf = lambda x : 0 if x < 0 else np.exp( -x )
>>> proposalCSampler = lambda x : np.random.normal( x, 1 )
>>> mhSampler = MetropolisHastingsSampler( initialVal, targetPdf,
...                                        proposalCsampler )

getSample()

Get a sample.

Returns: rst – Data point sample.
Return type: scalar or array_like of scalar

Examples

>>> sample = mhSampler.getSample()

Nataf transformation

This module implements the Nataf distribution. The Gaussian quadrature and root finding algorithm are used to determine the integral resutls and rho’ of Eq.(12) in reference 1. Reference 2 uses the Gauss–Legendre quadrature to calculate the integral results. With the transformation, Gauss–Hermite quadrature can also be used to calculate the integral results, as indicated in reference 3. However, the iteration method proposed in reference 3 requires Cholesky decomposition in the iteration. This module follows the Nataf transformation method implemented in reference 2.

1: Liu, P.L. and Der Kiureghian, A., 1986. Multivariate distribution models with prescribed marginals and covariances. Probabilistic engineering mechanics, 1(2), pp.105-112.
2(1,2): Ehre, M., Geyer, S., Kamariotis, A., Papaioannou, I., Sardi, L., 2022. Documentation of the ERA Distribution Classes.
3(1,2): Li, H., Lü, Z. and Yuan, X., 2008. Nataf transformation based point estimate method. Chinese Science Bulletin, 53(17), pp.2586-2592.

class ffpack.rpm.nataf.NatafTransformation(distObjs, corrMat, quadDeg=99, quadRange=8, randomSeed=None)

Nataf distribution for correlated marginal distributions.

__init__(distObjs, corrMat, quadDeg=99, quadRange=8, randomSeed=None)

Initialize the Nataf distribution.

Parameters

distObjs (array_like of distributions) – Marginal distribution objects. It should be the freezed distribution objects with pdf, cdf, ppf. We recommend to use scipy.stats functions.
corrMat (2d matrix) – Correlation matrix of the marginal distributions.
quadDeg (integer) – Quadrature degree.
quadRange (scalar) – Quadrature range. The integral will be performed in the range [ -quadRange, quadRange ].
randomSeed (integer, optional) – Random seed. If randomSeed is none or is not an integer, the random seed in global config will be used.

Raises

ValueError – If distObjs is empty. If dimensions are not match for distObjs and corrMat. If corrMat is not 2d matrix. If corrMat is not positive definite. If corrMat is not symmetric. If corrMat diagonal is not 1.

Examples

>>> from ffpack.rpm import NatafTransformation
>>> distObjs = [ stats.norm(), stats.norm() ]
>>> corrMat = [ [ 1.0, 0.5 ], [ 0.5, 1.0 ] ]
>>> natafDist = NatafTransformation( distObjs=distObjs, corrMat=corrMat )

cdf(X)

Get cdf value for X.

Parameters: X (1d array) – Data point in X space.
Returns: rst – Value of cumulative distribution function at data point X.
Return type: scalar
Raises: ValueError – If length of X does not match dim. If X is not 1d array.

Examples

>>> X = [ 0.5, 1.5 ]
>>> rst = natafDist.cdf( X )

getSample()

Get a sample in X space from Nataf distribution

Returns: X – Data point in X space.
Return type: 1d array

Examples

>>> X = natafDist.getSample()

getU(X)

Get data point in U space and Jacobian.

Parameters

X (1d array) – Data point in X space.

Returns

U (1d array) – Data point in U space.
J (2d matrix) – Jacobian from X space to U space.

Raises

ValueError – If length of X does not match dim. If X is not 1d array.

Notes

X -> Z -> U

Examples

>>> X = [ 0.5, 1.5 ]
>>> U, J = natafDist.getU( X )

getX(U)

Get data point in X space and Jacobian.

Parameters

U (1d array) – Data point in U space.

Returns

X (1d array) – Data point in X space.
J (2d matrix) – Jacobian from U space to X space.

Raises

ValueError – If length of U does not match dim. If U is not 1d array.

Notes

U -> Z -> X

Examples

>>> U = [ 0.5, 1.5 ]
>>> X, J = natafDist.getX( U )

pdf(X)

Get pdf value for X.

Parameters: X (1d array) – Data point in X space.
Returns: rst – Value of probability density function at data point X.
Return type: scalar
Raises: ValueError – If length of X does not match dim. If X is not 1d array.

Examples

>>> X = [ 0.5, 1.5 ]
>>> rst = natafDist.pdf( X )