RatioOfUniforms

(Source code, svg)

../../_images/RatioOfUniforms.svg
class RatioOfUniforms(*args)

Ratio of Uniforms simulation algorithm.

Warning

This class is experimental and likely to be modified in future releases. To use it, import the openturns.experimental submodule.

Available constructors:

RatioOfUniforms(dist)

RatioOfUniforms(logUnscaledF, rangeF, isScaled)

Parameters:
distDistribution

Distribution to sample.

With that constructor, the range is the distribution range and isScaled = True.

logUnscaledPDFFunction

Function \log F such that f = F\, 1_{I} is proportional to the PDF p we want to sample.

rangeFInterval

Interval I \subset \Rset^\inputDim outside which p=0.

isScaledbool, optional

Flag telling if the input logUnscaledPDF is such that f = F\, 1_{I} is an actual PDF of a distribution.

Default value is False.

Methods

asComposedEvent()

If the random vector can be viewed as the composition of several ThresholdEvent objects, this method builds and returns the composition.

getAcceptanceRatio()

Accessor to the acceptance ratio.

getAntecedent()

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

getC()

Accessor to the c normalization factor.

getCandidateNumber()

Accessor to the maximum number of points tested to find admissible starting points.

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.

getFrozenRealization(fixedPoint)

Compute realizations of the RandomVector.

getFrozenSample(fixedSample)

Compute realizations of the RandomVector.

getFunction()

Accessor to the Function in case of a composite RandomVector.

getInfV()

Accessor to the V-lower bound of \tilde{A}_{f,r}.

getLogUnscaledPDF()

Accessor to the function \log F.

getMarginal(*args)

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

getMaximumMultiStart()

Accessor to the maximum number of starting points for the multistart algorithm.

getMean()

Accessor to the mean of the RandomVector.

getName()

Accessor to the object's name.

getOperator()

Accessor to the comparaison operator of the Event.

getOptimizationAlgorithm()

Accessor to the optimization algorithm.

getParameter()

Accessor to the parameter of the distribution.

getParameterDescription()

Accessor to the parameter description of the distribution.

getProcess()

Get the stochastic process.

getR()

Accessor to the r parameter.

getRange()

Accessor to the interval I where the function f\propto p.

getRealization()

Generate a realization.

getSample(size)

Generate a sample.

getSampleWithTryNumber(size)

Sample with number of realizations needed to get it.

getSupU()

Accessor to the U-upper bound of \tilde{A}_{f,r}.

getSupV()

Accessor to the V-upper bound of \tilde{A}_{f,r}.

getThreshold()

Accessor to the threshold of the Event.

hasName()

Test if the object is named.

initialize()

Trigger to the initializationof the class.

isComposite()

Accessor to know if the RandomVector is a composite one.

isEvent()

Whether the random vector is an event.

isInitialized()

Flag telling if the algorithm has been initialized.

setCandidateNumber(candidateNumber)

Accessor to the maximum number of points tested to find admissible starting points.

setDescription(description)

Accessor to the description of the RandomVector.

setLogUnscaledPDFAndRange(logUnscaledPDF, ...)

Accessor to the function \log F and the range I.

setMaximumMultiStart(maximumMultiStart)

Accessor to the maximum number of starting points for the multistart algorithm.

setName(name)

Accessor to the object's name.

setOptimizationAlgorithm(optimizationAlgorithm)

Accessor to the optimization algorithm.

setParameter(parameters)

Accessor to the parameter of the distribution.

setR(r)

Accessor to the r parameter.

Notes

This algorithm samples a multivariate distribution whose probability density probability is known up to a multiplicative factor, see [stadlober1990] for the unidimensional case and [gobet2016] for the multivariate case.

Let p be a probability density function of dimension \inputDim and f: \Rset^\inputDim \rightarrow \Rset^+ a function such that:

p(\vect{x}) = cf(\vect{x})\quad  \forall \vect{x} \in \Rset^\inputDim

with c \in \Rset^+_*. We can write f as:

f(\vect{x}) = F(\vect{x})\, 1_{I}\quad \forall \vect{x}  \in I

with I \subset \Rset^\inputDim. It implies that the numerical range of p is inluced in I. We note:

I = \prod_{i=1}^\inputDim \left[a_i, b_i \right]

If the function f is an actual PDF and I its numerical range, then c=1. In the general case, c is unknown and assumed to be difficult to compute.

Let r \in \Rset^+_*. By default, r=1. Let A_{f,r} be the domain defined by:

A_{f,r} & = \left\{ (u, v_1, \dots, v_\inputDim) \in \Rset^{\inputDim+1}\, |\, 0 \leq u \leq f\left(\dfrac{v_1}{u^r},
\dots, \dfrac{v_\inputDim}{u^r}\right)^{\frac{1}{1+r\inputDim}}\right \}

The Lebesgue measure of A_{f,r} is equal to \dfrac{1}{c(1+rd)}.

Let (U, V_1, \dots, V_\inputDim) be a random variable uniformly distributed on A_{f,r}. Then, \left( \dfrac{V_1}{U^r}, \dots, \dfrac{V_\inputDim}{U^r}\right) is a random vector distributed according to p.

Under the condition that the functions \vect{x} \rightarrow f(\vect{x})^{\frac{1}{1+r\inputDim}} and \vect{x} \rightarrow x_if(\vect{x})^{\frac{1}{1+r\inputDim}} for all i are bounded, then the domain A_{f,r} is bounded and can be included in a rectangular bounding box \tilde{A}_{f,r} defined by:

(1)\tilde{A}_{f,r} & = \left[0, \sup_{\vect{x}} f(\vect{x})^{\frac{1}{1+r\inputDim}} \right] \times \prod_i
\left[ \inf_{\vect{x}} x_i f(\vect{x})^{\frac{r}{1+r\inputDim}} ,  \sup_{\vect{x}} x_i f(\vect{x})^{\frac{r}
{1+r\inputDim}}\right]\\
                & = [0, u_{sup}] \times \prod_{i=1}^d [v_{inf, i}, v_{sup, i}]

This allows one to sample uniformly the set A_{f,r} by rejection sampling inside \tilde{A}_{f,r}.

How to calculate the rectangular bounding box \tilde{A}_{f,r} ?

The upper bound u_{sup} is calculated by maximizing the log of the function rather than the function:

\sup_{\vect{x}} \log f(\vect{x})

To compute the second part of \tilde{A}_{f,r}, we note that:

  • If the upper bound b_i \leq 0, then:

    v_{sup, i} = \sup_{\vect{x}} x_i f(\vect{x})^{\frac{r}{1+rd}} =\sup_{\vect{x}, x_i \leq 0} x_i
f(\vect{x})^{\frac{r}{1+r\inputDim}} = 0

    Otherwise, we solve the maximization problem for x_i \in \Rset^+:

    v_{inf, i} = \sup_{\vect{x}} x_i f(\vect{x})^{\frac{r}{1+r\inputDim}} =\sup_{\vect{x}, x_i \geq 0} x_i
f(\vect{x})^{\frac{r}{1+r\inputDim}}

    which is equivalent to the problem:

    \sup_{\vect{x}, x_i \geq 0} \log x_i  + \dfrac{r}{1+r\inputDim} \log f(\vect{x})

  • If the lower bound a_i \geq 0, then:

    \inf_{\vect{x}} x_i f(\vect{x})^{\frac{r}{1+rd}} = \inf_{\vect{x}, x_i \geq 0} x_i
f(\vect{x})^{\frac{r}{1+r\inputDim}} = 0

    Otherwise, we solve the minimization problem for x_i \in \Rset^- and note that:

    \inf_{\vect{x} } x_i f(\vect{x})^{\frac{r}{1+r\inputDim}} = \inf_{\vect{x},  x_i \leq 0 } x_i
f(\vect{x})^{\frac{r}{1+r\inputDim}} = - \sup_{\vect{x},  x_i \leq 0 } |x_i| f(\vect{x})^{\frac{r
}{1+r\inputDim}}

    which is equivalent to the problem:

    \sup_{\vect{x}, x_i \leq 0} \log |x_i|  + \dfrac{r}{1+r\inputDim} \log f(\vect{x})

For the effective solution of these optimization problems, the initial point \vect{x}_0 must be such that f(\vect{x}_0)>0. Such points are found by scrolling through a Sobol multivariate low discrepancy sequence: the maximum number of points of the sequence to generate is fixed by the entry RatioOfUniforms-CandidateNumber of ResourceMap. The maximum number of points selected as starting points is controlled by the entry RatioOfUniforms-MaximumMultiStart. All the optimization problems are solved using a MultiStart algorithm based on these selected starting points. When no point has been selected, an exception is sent.

The acceptance ratio of the rejection algorithm is defined by:

(2)\tau_{accept} =  \frac{|A_{f,r}|}{|\tilde{A}_{f,r}|} = \dfrac{1}{c(1+r\inputDim)} \dfrac{1}{u_{sup}
\prod_{i=1}^d \left(v_{sup, i} - v_{inf, i} \right)}

It is computed by the initialize method: if c=1, the rate is computed by (2). If not, it is estimated as (3).

This class uses the following entries of ResourceMap:

  • RatioOfUniforms-OptimizationAlgorithm: the optimization algorithm used to compute the bounding box of \tilde{A}_{f,r},

  • RatioOfUniforms-CandidateNumber: maximum number of points of the Sobol low discrepancy sequence of dimension \inputDim used to find admissible starting points,

  • RatioOfUniforms-MaximumMultiStart: maximum number of points selected as starting points for the multistart optimization algorithm,

  • RatioOfUniforms-NormalizationSampleSize: size of the sample used to estimate the acceptance ratio.

Examples

Create a ratio of uniforms algorithm and get the multistart optimization problems used to compute the bounds of \tilde{A}_{f,r}:

>>> import openturns as ot
>>> import openturns.experimental as otexp
>>> from math import pi
>>> f = ot.SymbolicFunction('x', '(1.5+sin(x))*exp(x)')
>>> log_UnscaledPDF = ot.ComposedFunction(ot.SymbolicFunction('x', 'log(x)'), f)
>>> range_PDF = ot.Interval(0.0, 2.0 * pi)
>>> ratioOfU = otexp.RatioOfUniforms(log_UnscaledPDF, range_PDF, False)
>>> collMultiStart = ratioOfU.initialize()
>>> x = ratioOfU.getRealization()
>>> sample = ratioOfU.getSample(10)
>>> ratioOfU = otexp.RatioOfUniforms(ot.Student(8.5, 3))
__init__(*args)
asComposedEvent()

If the random vector can be viewed as the composition of several ThresholdEvent objects, this method builds and returns the composition. Otherwise throws.

Returns:
composedRandomVector

Composed event.

getAcceptanceRatio()

Accessor to the acceptance ratio.

Returns:
accepRatiofloat, 0 \leq \tau_{accept} \leq 1,

Acceptance ratio \tau_{accept}.

Notes

The acceptance ratio is computed at the initialization of the class by the ethod initialize, according to (2) where c=1 if the class has been created with a log-PDF and estimated by (3) in the other case.

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}).

getC()

Accessor to the c normalization factor.

Returns:
cfloat, c >0,

The c normalization factor.

Notes

If the class has been created with a scaled function, then c=1.

If the class has been created with an unscaled function, the acceptance ratio is estimated at the initialization of the class by the method initialize. From \hat{\tau}_{accept} defined in (3), we get an estimate of the normalization factor defined by:

\hat{c} = \dfrac{1}{\hat{\tau}_{accept}(1+r\inputDim)} \dfrac{1}{u_{sup}\prod_{i=1}^d \left(v_{sup, i} - v_{inf, i} \right)}

getCandidateNumber()

Accessor to the maximum number of points tested to find admissible starting points.

Returns:
nint, n>0

Maximum number of points tested to find admissible points for the optimization algorithm.

Notes

As long as the number of selected starting points is less than the number defined in ResourceMap, entry RatioOfUniforms-MaximumMultiStart, a new point of the Sobol low discrepancy sequence of dimension \inputDim is considered. The maximum points of the Sobol sequence is n. If no starting point has been selected, a warning is sent.

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.

getFrozenRealization(fixedPoint)

Compute realizations of the RandomVector.

In the case of a CompositeRandomVector or an event of some kind, this method returns the value taken by the random vector if the root cause takes the value given as argument.

Parameters:
fixedPointPoint

Point chosen as the root cause of the random vector.

Returns:
realizationPoint

The realization corresponding to the chosen root cause.

Examples

>>> import openturns as ot
>>> distribution = ot.Normal()
>>> randomVector = ot.RandomVector(distribution)
>>> f = ot.SymbolicFunction('x', 'x')
>>> compositeRandomVector = ot.CompositeRandomVector(f, randomVector)
>>> event = ot.ThresholdEvent(compositeRandomVector, ot.Less(), 0.0)
>>> print(event.getFrozenRealization([0.2]))
[0]
>>> print(event.getFrozenRealization([-0.1]))
[1]
getFrozenSample(fixedSample)

Compute realizations of the RandomVector.

In the case of a CompositeRandomVector or an event of some kind, this method returns the different values taken by the random vector when the root cause takes the values given as argument.

Parameters:
fixedSampleSample

Sample of root causes of the random vector.

Returns:
sampleSample

Sample of the realizations corresponding to the chosen root causes.

Examples

>>> import openturns as ot
>>> distribution = ot.Normal()
>>> randomVector = ot.RandomVector(distribution)
>>> f = ot.SymbolicFunction('x', 'x')
>>> compositeRandomVector = ot.CompositeRandomVector(f, randomVector)
>>> event = ot.ThresholdEvent(compositeRandomVector, ot.Less(), 0.0)
>>> print(event.getFrozenSample([[0.2], [-0.1]]))
    [ y0 ]
0 : [ 0  ]
1 : [ 1  ]
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}).

getInfV()

Accessor to the V-lower bound of \tilde{A}_{f,r}.

Returns:
infVPoint

The V-lower bound of \tilde{A}_{f,r} defined in (1).

Notes

The V-lower bound of \tilde{A}_{f,r} is \left( \inf_{\vect{x}} x_i f(\vect{x})^{\frac{r}{1+r\inputDim}} \right)_{1 \leq i \leq \inputDim}.

getLogUnscaledPDF()

Accessor to the function \log F.

Returns:
logUnscaledPDFFunction,

Function \log F such that f = F\, 1_{I} is proportional to the PDF p we want to sample.

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)
getMaximumMultiStart()

Accessor to the maximum number of starting points for the multistart algorithm.

Returns:
nint

Maximum number of starting points for the multistart algorithm.

Notes

A tested point \vect{x} \in I is a starting point of the multistart algorithm if f(\vect{x})>0.

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]
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.

getOptimizationAlgorithm()

Accessor to the optimization algorithm.

Returns:
algoOptimizationAlgorithm

Optimization algorithm that computes the bounds of \tilde{A}_{f,r}.

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.

getR()

Accessor to the r parameter.

Returns:
rfloat

The r parameter.

getRange()

Accessor to the interval I where the function f\propto p.

Returns:
rangeFInterval

Interval I \subset \Rset^\inputDim where the function f\propto p.

Notes

As p \propto f = F\, 1_I, then p = 0 outside I.

getRealization()

Generate a realization.

Returns:
point: Point

A pseudo-random realization of the distribution.

getSample(size)

Generate a sample.

Returns:
sample: Sample

A pseudo-random sample of the distribution.

getSampleWithTryNumber(size)

Sample with number of realizations needed to get it.

Parameters:
sizeint,

The number of accepted realizations of the rejection algorithm.

Returns:
sample: Sample,

A pseudo-random realization of the distribution of size size.

tryNumberint,

The number realizations needed to produce size accepted realizations.

Notes

The rejection algorithm produces tryNumber realizations of the uniform distribution on \tilde{A}_{f,r} to produce size realizations inside A_{f,r}.

getSupU()

Accessor to the U-upper bound of \tilde{A}_{f,r}.

Returns:
supUfloat, positive

The U-upper bound of \tilde{A}_{f,r} defined in (1).

Notes

The U-upper bound of \tilde{A}_{f,r} is \sup_{\vect{x}} f(\vect{x})^{\frac{1}{1+r\inputDim}}.

getSupV()

Accessor to the V-upper bound of \tilde{A}_{f,r}.

Returns:
infVPoint

The V-upper bound of \tilde{A}_{f,r} defined in (1).

Notes

The V-upper bound of \tilde{A}_{f,r} is \left( \sup_{\vect{x}} x_i f(\vect{x})^{\frac{r}{1+r\inputDim}} \right)_{1 \leq i \leq \inputDim}.

getThreshold()

Accessor to the threshold of the Event.

Returns:
thresholdfloat

Threshold of the RandomVector.

hasName()

Test if the object is named.

Returns:
hasNamebool

True if the name is not empty.

initialize()

Trigger to the initializationof the class.

Returns:
collMutliStartcollection of MultiStart

Collection of the multistart optimization algorithms used to compute the bounds of \tilde{A}_{f,r}.

Notes

A call to that method triggers the computation of the bounding box \tilde{A}_{f,r} and the estimation of the acceptance ratio.

The methods returns a collection of multistart optimization algorithms containing:

  • the multistart optimization algorithm used to compute u_{sup},

  • for each i, the multistart optimization algorithm used to compute v_{inf, i} if required, the one used to compute v_{sup, i} if required.

Then, the size of the collection of MultiStart is less than 2\inputDim + 1 and at least 1. These multistart optimization algorithms give access to the starting points used to compute the optimization problem and all the results associated.

Within the initialization of the class, the acceptance ratio is computed:

  • If the class has been created with a scaled function, the acceptance ratio is defined by (2) where c=1. If the acceptance ratio is greater than 1, then it means that the optimization problem has not worked well and a warning is sent.

  • If the class has been created with an unscaled function, the acceptance ratio is estimated with the method getSampleWithTryNumber: if n_{try} is the number of realizations uniformly generated inside \tilde{A}_{f,r} to get n_{accept} accepted realizations, then we have:

    (3)\hat{\tau}_{accept} = \dfrac{n_{accept}}{n_{try}}

    In that case, \hat{\tau}_{accept} is always less than 1.

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}.

isInitialized()

Flag telling if the algorithm has been initialized.

Returns:
answer: bool

Flag telling if the bounding box \tilde{A}_{f,r} and the acceptance ratio has already been computed.

setCandidateNumber(candidateNumber)

Accessor to the maximum number of points tested to find admissible starting points.

Parameters:
nint, n>0

Maximum number of points tested to find admissible points for the optimization algorithm.

Notes

As long as the number of selected starting points is less than the number defined in ResourceMap, entry RatioOfUniforms-MaximumMultiStart, a new point of the Sobol low discrepancy sequence of dimension \inputDim is considered. The maximum points of the Sobol sequence is n. If no starting point has been selected, a warning is sent.

setDescription(description)

Accessor to the description of the RandomVector.

Parameters:
descriptionstr or sequence of str

Describes the components of the RandomVector.

setLogUnscaledPDFAndRange(logUnscaledPDF, range, isScaled)

Accessor to the function \log F and the range I.

Parameters:
logUnscaledPDFFunction

Function \log F where F\, 1_{I} \propto p, the PDF we want to sample.

rangeFInterval

Interval I \in \Rset^\inputDim outside which p=0.

setMaximumMultiStart(maximumMultiStart)

Accessor to the maximum number of starting points for the multistart algorithm.

Parameters:
nint, n>0

Maximum number of starting points for the multistart algorithm.

Notes

A tested point \vect{x} \in I is kept as a starting point of the multistart algorithm if f(\vect{x})>0.

setName(name)

Accessor to the object’s name.

Parameters:
namestr

The name of the object.

setOptimizationAlgorithm(optimizationAlgorithm)

Accessor to the optimization algorithm.

Parameters:
algoOptimizationAlgorithm

Optimization algorithm that computes the bounds of \tilde{A}_{f,r}.

setParameter(parameters)

Accessor to the parameter of the distribution.

Parameters:
parametersequence of float

Parameter values.

setR(r)

Accessor to the r parameter.

Parameters:
rfloat

The r parameter.

Notes

This method triggers the initialization of the class.

Examples using the class

Sampling from an unscaled probability density

Sampling from an unscaled probability density

Use the Ratio of Uniforms algorithm to sample a distribution

Use the Ratio of Uniforms algorithm to sample a distribution