Gibbs¶

class Gibbs(*args)¶

Gibbs sampling method.

Refer to Bayesian calibration, The Metropolis-Hastings Algorithm.

Parameters:

samplerssequence of MetropolisHastings: List of samplers for each bloc of the chain

See also

RandomWalkMetropolisHastings, RandomVectorMetropolisHastings

Notes

Assume we want to sample from a complicated joint distribution. The Gibbs algorithm is a Markov Chain Monte-Carlo algorithm which calls samplers (implemented as MetropolisHastings objects) in a predetermined sequence (by default) or randomly: see setUpdatingMethod() to manually choose the behavior, and note the default behavior can be changed by setting the Gibbs-DefaultUpdatingMethod ResourceMap key. Each sampler samples from the conditional distributions of one or several components.

Examples

>>> import openturns as ot
>>> ot.RandomGenerator.SetSeed(0)
>>> chainDim = 3
>>> # Observations
>>> obsDim = 1
>>> obsSize = 10
>>> y = [-9.50794871493506, -3.83296694500105, -2.44545713047953,
...      0.0803625289211318, 1.01898069723583, 0.661725805623086,
...      -1.57581204592385, -2.95308465670895, -8.8878164296758,
...      -13.0812290405651]
>>> y_obs = ot.Sample([[yi] for yi in y])
>>> # Parameters
>>> covariates = ot.Sample(obsSize, chainDim)
>>> for i in range(obsSize):
...     for j in range(chainDim):
...         covariates[i, j] = (-2 + 5.0 * i / 9.0) ** j
>>> # Model
>>> fullModel = ot.SymbolicFunction(
...          ['p1', 'p2', 'p3', 'x1', 'x2', 'x3'],
...          ['p1*x1+p2*x2+p3*x3', '1.0'])
>>> parametersSet = range(chainDim)
>>> parametersValue = [0.0] * len(parametersSet) # 0.0 is a placeholder
>>> linkFunction = ot.ParametricFunction(fullModel, parametersSet, parametersValue)
>>> # Calibration parameters
>>> # Proposal distribution
>>> prop = ot.Uniform(-1.0, 1.0)
>>> # Prior distribution
>>> sigma0 = [10.0]*chainDim
>>> #  Covariance matrix
>>> Q0_inv = ot.CorrelationMatrix(chainDim)
>>> for i in range(chainDim): 
...     Q0_inv[i, i] = sigma0[i] * sigma0[i]
>>> mu0 = [0.0]*chainDim
>>> # x0 ~ N(mu0, sigma0)
>>> prior = ot.Normal(mu0, Q0_inv)
>>> # Conditional distribution y~N(z, 1.0)
>>> conditional = ot.Normal()
>>> # Create a Gibbs sampler
>>> coll = [ot.RandomWalkMetropolisHastings(prior, mu0, prop, [i]) for i in range(chainDim)]
>>> for mh in coll: mh.setLikelihood(conditional, y_obs, linkFunction, covariates)
>>> sampler = ot.Gibbs(coll)
>>> sampler.setBurnIn(20)
>>> sampler.setThinning(10)
>>> # Get a realization
>>> mu = sampler.getRealization()
>>> # Create a Gibbs sampler which updates one randomly chosen component at each step
>>> sampler2 = ot.Gibbs(coll)
>>> sampler2.setUpdatingMethod(ot.Gibbs.RANDOM_UPDATING)
>>> mu2 = sampler2.getRealization()

Methods

`getAntecedent`()	Accessor to the antecedent RandomVector in case of a composite RandomVector.
`getBurnIn`()	Get the length of the burn-in period.
`getClassName`()	Accessor to the object's name.
`getCovariance`()	Accessor to the covariance of the RandomVector.
`getDescription`()	Accessor to the description of the RandomVector.
`getDimension`()	Accessor to the dimension of the RandomVector.
`getDistribution`()	Accessor to the distribution of the RandomVector.
`getDomain`()	Accessor to the domain of the Event.
`getFunction`()	Accessor to the Function in case of a composite RandomVector.
`getHistory`()	Get the history storage.
`getId`()	Accessor to the object's id.
`getMarginal`(*args)	Get the random vector corresponding to the $i^{th}$ marginal component(s).
`getMean`()	Accessor to the mean of the RandomVector.
`getMetropolisHastingsCollection`()	Get the list of Metropolis-Hastings samplers.
`getName`()	Accessor to the object's name.
`getOperator`()	Accessor to the comparaison operator of the Event.
`getParameter`()	Accessor to the parameter of the distribution.
`getParameterDescription`()	Accessor to the parameter description of the distribution.
`getProcess`()	Get the stochastic process.
`getRealization`()	Compute one realization of the RandomVector.
`getRecomputeLogPosterior`()	Get the posterior recompute flags.
`getSample`(size)	Compute realizations of the RandomVector.
`getShadowedId`()	Accessor to the object's shadowed id.
`getThinning`()	Get the thinning parameter.
`getThreshold`()	Accessor to the threshold of the Event.
`getUpdatingMethod`()	Get how the order of the samplers is determined.
`getVisibility`()	Accessor to the object's visibility state.
`hasName`()	Test if the object is named.
`hasVisibleName`()	Test if the object has a distinguishable name.
`isComposite`()	Accessor to know if the RandomVector is a composite one.
`isEvent`()	Whether the random vector is an event.
`setBurnIn`(burnIn)	Set the length of the burn-in period.
`setDescription`(description)	Accessor to the description of the RandomVector.
`setHistory`(strategy)	Set the history storage.
`setName`(name)	Accessor to the object's name.
`setParameter`(parameters)	Accessor to the parameter of the distribution.
`setShadowedId`(id)	Accessor to the object's shadowed id.
`setThinning`(thinning)	Set the thinning parameter.
`setUpdatingMethod`(updatingMethod)	Set how the order of the samplers is determined.
`setVisibility`(visible)	Accessor to the object's visibility state.

__init__(*args)¶

getAntecedent()¶

Accessor to the antecedent RandomVector in case of a composite RandomVector.

Returns:

antecedentRandomVector: Antecedent RandomVector $\vect{X}$ in case of a CompositeRandomVector such as: $\vect{Y}=f(\vect{X})$ .

getBurnIn()¶

Get the length of the burn-in period.

Returns:

lengthint: Length of the burn-in period, that is the number of first iterates of the MCMC chain which will be thrown away when generating the sample.

getClassName()¶

Accessor to the object’s name.

Returns:

class_namestr: The object class name (object.__class__.__name__).

getCovariance()¶

Accessor to the covariance of the RandomVector.

Returns:

covarianceCovarianceMatrix: Covariance of the considered UsualRandomVector.

Examples

>>> import openturns as ot
>>> distribution = ot.Normal([0.0, 0.5], [1.0, 1.5], ot.CorrelationMatrix(2))
>>> randomVector = ot.RandomVector(distribution)
>>> ot.RandomGenerator.SetSeed(0)
>>> print(randomVector.getCovariance())
[[ 1    0    ]
 [ 0    2.25 ]]

getDescription()¶

Accessor to the description of the RandomVector.

Returns:

descriptionDescription: Describes the components of the RandomVector.

getDimension()¶

Accessor to the dimension of the RandomVector.

Returns:

dimensionpositive int: Dimension of the RandomVector.

getDistribution()¶

Accessor to the distribution of the RandomVector.

Returns:

distributionDistribution: Distribution of the considered UsualRandomVector.

Examples

>>> import openturns as ot
>>> distribution = ot.Normal([0.0, 0.0], [1.0, 1.0], ot.CorrelationMatrix(2))
>>> randomVector = ot.RandomVector(distribution)
>>> ot.RandomGenerator.SetSeed(0)
>>> print(randomVector.getDistribution())
Normal(mu = [0,0], sigma = [1,1], R = [[ 1 0 ]
 [ 0 1 ]])

getDomain()¶

Accessor to the domain of the Event.

Returns:

domainDomain: Describes the domain of an event.

getFunction()¶

Accessor to the Function in case of a composite RandomVector.

Returns:

functionFunction: Function used to define a CompositeRandomVector as the image through this function of the antecedent $\vect{X}$ : $\vect{Y}=f(\vect{X})$ .

getHistory()¶

Get the history storage.

Returns:

historyHistoryStrategy: Used to record the chain.

getId()¶

Accessor to the object’s id.

Returns:

idint: Internal unique identifier.

getMarginal(*args)¶

Get the random vector corresponding to the $i^{th}$ marginal component(s).

Parameters:

iint or list of ints, $0\leq i < dim$: Indicates the component(s) concerned. $dim$ is the dimension of the RandomVector.

Returns:

vectorRandomVector: RandomVector restricted to the concerned components.

Notes

Let’s note $\vect{Y}=\Tr{(Y_1,\dots,Y_n)}$ a random vector and $I \in [1,n]$ a set of indices. If $\vect{Y}$ is a UsualRandomVector, the subvector is defined by $\tilde{\vect{Y}}=\Tr{(Y_i)}_{i \in I}$ . If $\vect{Y}$ is a CompositeRandomVector, defined by $\vect{Y}=f(\vect{X})$ with $f=(f_1,\dots,f_n)$ , $f_i$ some scalar functions, the subvector is $\tilde{\vect{Y}}=(f_i(\vect{X}))_{i \in I}$ .

Examples

>>> import openturns as ot
>>> distribution = ot.Normal([0.0, 0.0], [1.0, 1.0], ot.CorrelationMatrix(2))
>>> randomVector = ot.RandomVector(distribution)
>>> ot.RandomGenerator.SetSeed(0)
>>> print(randomVector.getMarginal(1).getRealization())
[0.608202]
>>> print(randomVector.getMarginal(1).getDistribution())
Normal(mu = 0, sigma = 1)

getMean()¶

Accessor to the mean of the RandomVector.

Returns:

meanPoint: Mean of the considered UsualRandomVector.

Examples

>>> import openturns as ot
>>> distribution = ot.Normal([0.0, 0.5], [1.0, 1.5], ot.CorrelationMatrix(2))
>>> randomVector = ot.RandomVector(distribution)
>>> ot.RandomGenerator.SetSeed(0)
>>> print(randomVector.getMean())
[0,0.5]

getMetropolisHastingsCollection()¶

Get the list of Metropolis-Hastings samplers.

Returns:

samplerssequence of MetropolisHastings: List of all Metropolis-Hastings samplers used in the Gibbs algorithm

getName()¶

Accessor to the object’s name.

Returns:

namestr: The name of the object.

getOperator()¶

Accessor to the comparaison operator of the Event.

Returns:

operatorComparisonOperator: Comparaison operator used to define the RandomVector.

getParameter()¶

Accessor to the parameter of the distribution.

Returns:

parameterPoint: Parameter values.

getParameterDescription()¶

Accessor to the parameter description of the distribution.

Returns:

descriptionDescription: Parameter names.

getProcess()¶

Get the stochastic process.

Returns:

processProcess: Stochastic process used to define the RandomVector.

getRealization()¶

Compute one realization of the RandomVector.

Returns:

aRealizationPoint: Sequence of values randomly determined from the RandomVector definition. In the case of an event: one realization of the event (considered as a Bernoulli variable) which is a boolean value (1 for the realization of the event and 0 else).

See also

getRealization

Examples

>>> import openturns as ot
>>> distribution = ot.Normal([0.0, 0.0], [1.0, 1.0], ot.CorrelationMatrix(2))
>>> randomVector = ot.RandomVector(distribution)
>>> ot.RandomGenerator.SetSeed(0)
>>> print(randomVector.getSample(3))
    [ X0        X1        ]
0 : [  0.608202 -1.26617  ]
1 : [ -0.438266  1.20548  ]
2 : [ -2.18139   0.350042 ]

getShadowedId()¶

Accessor to the object’s shadowed id.

Returns:

idint: Internal unique identifier.

getThinning()¶

Get the thinning parameter.

Returns:

thinningint: Thinning parameter: storing only every $k^{th}$ point after the burn-in period.

Notes

When generating a sample of size $q$ , the number of MCMC iterations performed is $l+1+(q-1)k$ where $l$ is the burn-in period length and $k$ the thinning parameter.

getThreshold()¶

Accessor to the threshold of the Event.

Returns:

thresholdfloat: Threshold of the RandomVector.

getUpdatingMethod()¶

Get how the order of the samplers is determined.

Available values are:

0 (Gibbs.DETERMINISTIC_UPDATING): the samplers are sequentially called in the given order.

1 (Gibbs.RANDOM_UPDATING): one of the samplers is randomly selected to produce the next realization.

Refer to the setUpdatingMethod() documentation for details.

Returns:

updatingMethodint: See above for the possible values and their meaning.

getVisibility()¶

Accessor to the object’s visibility state.

Returns:

visiblebool: Visibility flag.

hasName()¶

Test if the object is named.

Returns:

hasNamebool: True if the name is not empty.

hasVisibleName()¶

Test if the object has a distinguishable name.

Returns:

hasVisibleNamebool: True if the name is not empty and not the default one.

isComposite()¶

Accessor to know if the RandomVector is a composite one.

Returns:

isCompositebool: Indicates if the RandomVector is of type Composite or not.

isEvent()¶

Whether the random vector is an event.

Returns:

isEventbool: Whether it takes it values in {0, 1}.

setBurnIn(burnIn)¶

Set the length of the burn-in period.

Parameters:

lengthint: Length of the burn-in period, that is the number of first iterates of the MCMC chain which will be thrown away when generating the sample.

setDescription(description)¶

Accessor to the description of the RandomVector.

Parameters:

descriptionstr or sequence of str: Describes the components of the RandomVector.

setHistory(strategy)¶

Set the history storage.

Parameters:

historyHistoryStrategy: Used to record the chain.

setName(name)¶

Accessor to the object’s name.

Parameters:

namestr: The name of the object.

setParameter(parameters)¶

Accessor to the parameter of the distribution.

Parameters:

parametersequence of float: Parameter values.

setShadowedId(id)¶

Accessor to the object’s shadowed id.

Parameters:

idint: Internal unique identifier.

setThinning(thinning)¶

Set the thinning parameter.

Parameters:

thinningint, $k \geq 0$: Thinning parameter: storing only every $k^{th}$ point after the burn-in period.

Notes

When generating a sample of size $q$ , the number of MCMC iterations performed is $l+1+(q-1)k$ where $l$ is the burn-in period length and $k$ the thinning parameter.

setUpdatingMethod(updatingMethod)¶

Set how the order of the samplers is determined.

Available values are:

0 (Gibbs.DETERMINISTIC_UPDATING): the samplers are sequentially called in the given order.

1 (Gibbs.RANDOM_UPDATING): one of the samplers is randomly selected to produce the next realization.

Note that under DETERMINISTIC_UPDATING, all samplers are called in the specified order every time getRealization() is called. By contrast, under RANDOM_UPDATING, only one sampler is called by getRealization(), and this sampler is randomly chosen among all specified samplers.

Parameters: