Note

Go to the end to download the full example code.

Use the Ratio of Uniforms algorithm to sample a distribution¶

In this example, we illustrate in dimension $\inputDim = 1$ the ratio of uniforms algorithm to sample a distribution which probability density function is known up to a normalization factor.

Consider a scalar distribution whose probability density function $p$ is known up to a normalization factor: let $f: \Rset \rightarrow \Rset^+$ be a function such that $p(x) = cf(x)$ for all $x \in \Rset$ with $c \in \Rset^+_*$ .

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

$A_{f,r} = \left\{ (u, v) \in \Rset^2 \, |\, 0 \leq u \leq f\left(\dfrac{v}{u^r}\right)^{\frac{1}{1+r}}\right \}$

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

Let $(U, V)$ be a random variable uniformly distributed on $A_{f,r}$ . Then, $X = \dfrac{V}{U^r}$ is a random variable on $\Rset$ distributed according to $p$ .

Under the condition that $f(x)^{\frac{1}{1+r}}$ and $x f(x)^{\frac{1}{1+r}}$ 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:

$\tilde{A}_{f,r} & = \left[0, \sup_{x} f(x)^{\frac{1}{1+r}} \right] \times \left[ \inf_{x} x f(x)^{\frac{1}{1+r}} , \sup_{x} x f(x)^{\frac{1}{1+r}}\right]\\ & = \left[0, u_{sup} \right] \times \left[v_{inf}, v_{sup}\right]$

Note that $u_{sup} > 0$ , $v_{sup} \geq 0$ and $v_{inf} \leq 0$ .

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

Study of the $(u,v)$ -space.¶

The boundary of $\tilde{A}_{f,r}$ is defined by:

$\partial \tilde{A}_{f,r} = \cF_1 \cup \cF_2 \cup \cF_3 \cup \cF_4$

where the boundaries $\cF_i$ are defined by:

$\cF_1 & = \left \{ (u, v_{sup}), 0 \leq u \leq u_{sup} \right \} \\ \cF_2 & = \left \{ (u, v_{inf}), 0 \leq u \leq u_{sup} \right \}\\ \cF_3 & = \left \{ (u_{sup}, v), v_{inf} \leq v \leq v_{sup} \right \}\\ \cF_4 & = \left \{ (0, v), v_{inf} \leq v \leq v_{sup} \right \}$

The boundary of $A_{f,r}$ is defined by:

$\partial A_{f,r} = \left\{ (u,v) = \left( f(x)^{\dfrac{1}{1+r}}, x f(x)^{\dfrac{r}{1+r}} \right), x \in \Rset \right\}$

The following function draws the domain $A_{f,r}$ and the bounding box $\partial \tilde{A}_{f,r}$ in the $(u,v)$ -space.

import openturns as ot
import openturns.viewer as otv
import openturns.experimental as otexp
from math import pi, exp


def draw_RatioOfUniforms_UV(ratio_algo):
    r = ratio_algo.getR()
    range_logF = ratio_algo.getRange()
    xInf = range_logF.getLowerBound()[0]
    xSup = range_logF.getUpperBound()[0]
    uSup = ratio_algo.getSupU()
    vInf = ratio_algo.getInfV()[0]
    vSup = ratio_algo.getSupV()[0]
    logF = ratio_algo.getLogUnscaledPDF()

    # \partial A_{f,r} domain
    x = xInf
    u = 0.0
    v = 0.0
    xPrec = None
    uPrec = None
    vPrec = None
    first = True
    step = 0.01 * (uSup + vSup - vInf)
    delta = step
    data_UV = ot.Sample(0, 2)
    while x < xSup:
        val_F = exp(logF([x])[0])
        u = val_F ** (1.0 / (1.0 + r))
        v = x * val_F ** (r / (1.0 + r))
        # Step adaptation after the first point
        if not first:
            # If the step is too large
            while abs(u - uPrec) + abs(v - vPrec) > 2 * delta:
                step *= 0.5
                x = xPrec + step
                if x < xSup:
                    val_F = exp(logF([x])[0])
                    u = val_F ** (1.0 / (1.0 + r))
                    v = x * val_F ** (r / (1.0 + r))
            # If the step is too small
            while abs(u - uPrec) + abs(v - vPrec) < 0.5 * delta and x < xSup:
                step *= 2.0
                x = xPrec + step
                if x < xSup:
                    val_F = exp(logF([x])[0])
                    u = val_F ** (1.0 / (1.0 + r))
                    v = x * val_F ** (r / (1.0 + r))
        xPrec = x
        uPrec = u
        vPrec = v
        x += step
        first = False
        data_UV.add([u, v])
    # to get a closed boundary
    data_UV.add(data_UV[0])
    g = ot.Graph("", "u", "v", True, "bottomright")
    data_A = ot.Curve(data_UV)
    data_A.setColor("purple")
    g.add(data_A)

    # \partial \tilde{A}_{f,r} bounding box
    boundary_1 = ot.Curve([[0, vSup], [uSup, vSup]])
    boundary_1.setColor("blue")
    boundary_2 = ot.Curve([[0, vInf], [uSup, vInf]])
    boundary_2.setColor("red")
    boundary_3 = ot.Curve([[uSup, vInf], [uSup, vSup]])
    boundary_3.setColor("green")
    boundary_4 = ot.Curve([[0, vInf], [0, vSup]])
    boundary_4.setColor("black")
    g.add(boundary_1)
    g.add(boundary_2)
    g.add(boundary_3)
    g.add(boundary_4)

    leg = g.getLegends()
    leg[0] = r"boundary of $A_{f,r}$"
    leg[1] = r"$\mathcal{F}_1$"
    leg[2] = r"$\mathcal{F}_2$"
    leg[3] = r"$\mathcal{F}_3$"
    leg[4] = r"$\mathcal{F}_4$"
    g.setLegends(leg)

    return g

Study of the $(x,y)$ -space.¶

We use the mapping function $\varphi$ defined by:

$\varphi: & \Rset^2 \rightarrow \Rset^2 \\ & (u,v) \rightarrow \left( \dfrac{v}{u^r}, u \right)$

Then we have:

$\varphi(A_{f,r}) = \left\{ (x,y) \, |\, x \in \Rset, 0 \leq y \leq f(x)^{\dfrac{1}{1+r}} \right\}$

and:

$\partial \varphi(\tilde{A}_{f,r}) = \varphi(\cF_1) \cup \varphi(\cF_2) \cup \varphi(\cF_3)$

where the boundaries $\varphi(\cF_i)$ are defined by:

$\varphi(\cF_1) & = \left\{ \varphi(u, v_{sup}), 0 \leq u \leq u_{sup} \right\}\\ & = \left\{ \left(x,\left(\dfrac{v_{sup}}{x}\right)^{1/r} \right), x \geq \dfrac{v_{sup}}{u_{sup}^{r}} \right \}$

$\varphi(\cF_2) & = \left \{ \varphi(u, v_{inf}), 0 \leq u \leq u_{sup} \right\}\\ & = \left\{ \left(x,\left(\dfrac{v_{inf}}{x}\right)^{1/r}\right), x \leq \dfrac{v_{inf}}{u_{sup}^{r}} \right \}$

$\varphi(\cF_3) & = \left\{ \varphi(u_{sup}, v), v_{inf} \leq v \leq v_{sup} \right\}\\ & = \left\{ (x, u_{sup}), \dfrac{v_{inf}}{u_{sup}^{r}} \leq x \leq \dfrac{v_{sup}}{u_{sup}^{r}} \right\}$

The following function draws the boundary $\partial \varphi(A_{f,r})$ and the bounding box $\partial \tilde{A}_{f,r}$ in the $(x,y)$ -space.

def draw_RatioOfUniforms_XY(ratio_algo, nPoints):
    r = ratio_algo.getR()
    range_logF = ratio_algo.getRange()
    xInf = range_logF.getLowerBound()[0]
    xSup = range_logF.getUpperBound()[0]
    uSup = ratio_algo.getSupU()
    vInf = ratio_algo.getInfV()[0]
    vSup = ratio_algo.getSupV()[0]
    logF = ratio_algo.getLogUnscaledPDF()

    # \partial varpĥi(A_{f,r}) domain
    eps = 1e-10
    step = (xSup - xInf - 2.0 * eps) / (nPoints - 1)
    data_XY = ot.Sample(nPoints, 2)
    for i in range(nPoints):
        x = xInf + eps + i * step
        val_F = exp(logF([x])[0])
        data_XY[i, 0] = x
        data_XY[i, 1] = val_F ** (1.0 / (1.0 + r))
    g = ot.Graph("", "x", "y", True, "lower center")
    data_A = ot.Curve(data_XY)
    data_A.setColor("purple")
    g.add(data_A)

    # \partial varphi(\tilde{A}_{f,r}) bounding box
    leg = ot.Description()
    leg.add(r"boundary of $\varphi(A_{f,r})$")
    boundary_1_fct = ot.SymbolicFunction(
        "x", "(" + str(vSup) + "/x)^(1/" + str(r) + ")"
    )
    if vSup / (uSup**r) < xSup:
        boundary_1 = boundary_1_fct.draw(vSup / (uSup**r), xSup, 101)
        boundary_1.setColors(["blue"])
        leg.add(r"$\varphi(\mathcal{F}_1)$")
        g.add(boundary_1)
    if vInf / (uSup**r) < uSup:
        boundary_2_fct = ot.SymbolicFunction(
            "x", "(" + str(vInf) + "/x)^(1/" + str(r) + ")"
        )
        boundary_2 = boundary_2_fct.draw(xInf, vInf / (uSup**r), 101)
        boundary_2.setColors(["red"])
        leg.add(r"$\varphi(\mathcal{F}_2)$")
        g.add(boundary_2)
    boundary_3 = ot.Curve([[vInf / (uSup**r), uSup], [vSup / (uSup**r), uSup]])
    boundary_3.setColor("green")
    g.add(boundary_3)
    leg.add(r"$\varphi(\mathcal{F}_3)$")
    g.setLegends(leg)

    return g

Illustration.¶

We create the function $\log f$ defined by:

(1)¶ $\log f(x) = \log(\cos(x)) + x, \quad x \in \left[-\dfrac{\pi}{2}, \dfrac{\pi}{2}\right]$

and the RatioOfUniforms.

log_UnscaledPDF = ot.SymbolicFunction("x", "log(cos(x)) + x")
eps = 1e-5
range_PDF = ot.Interval(-pi / 2.0 + eps, pi / 2.0 - eps)
ratio_Algo = otexp.RatioOfUniforms(log_UnscaledPDF, range_PDF)

We draw the boundary $\partial A_{f,r}$ and the bounding box $\partial \tilde{A}_{f,r}$ in the $(u,v)$ -space.

g_UV = draw_RatioOfUniforms_UV(ratio_Algo)
g_UV.setTitle(
    rf"$r = ${ratio_Algo.getR():.3g}, $\tau = ${ratio_Algo.getAcceptanceRatio():.3g}"
)

We draw the boundary $\partial \varphi(A_{f,r})$ and the bounding box $\partial \varphi(\tilde{A}_{f,r})$ in the $(x,y)$ -space. sphinx_gallery_thumbnail_number = 1

g_XY = draw_RatioOfUniforms_XY(ratio_Algo, 1001)
g_XY.setTitle(
    rf"$r = ${ratio_Algo.getR():.3g}, $\tau = ${ratio_Algo.getAcceptanceRatio():.3g}"
)

grid_default = ot.GridLayout(1, 2)
grid_default.setGraph(0, 0, g_UV)
grid_default.setGraph(0, 1, g_XY)
grid_default.setTitle(r"$\log f(x) = \log(\cos(x)) + x, x \in ]-\pi/2, \pi/2[$")
view = otv.View(grid_default)

$$\log f(x) = \log(\cos(x)) + x, x \in ]-\pi/2, \pi/2[$, $r = $1, $\tau = $0.564, $r = $1, $\tau = $0.564$

We change the $r$ value and we draw the same plots.

grid_full = ot.GridLayout(4, 2)
for i, r in enumerate([0.5, 0.8, 1.2, 2.0]):
    ratio_Algo.setR(r)
    g_UV = draw_RatioOfUniforms_UV(ratio_Algo)
    g_UV.setLegendPosition("")
    g_UV.setTitle(
        rf"$r = ${ratio_Algo.getR():.3g}, $\tau = ${ratio_Algo.getAcceptanceRatio():.3g}"
    )

    g_XY = draw_RatioOfUniforms_XY(ratio_Algo, 1001)
    g_XY.setLegendPosition("")
    g_XY.setTitle(
        rf"$r = ${ratio_Algo.getR():.3g}, $\tau = ${ratio_Algo.getAcceptanceRatio():.3g}"
    )

    grid_full.setGraph(i, 0, g_UV)
    grid_full.setGraph(i, 1, g_XY)

grid_full.setTitle(r"$\log f(x) = \log(\cos(x)) + x,  \quad x \in ]-\pi/2, \pi/2[$")
view = otv.View(grid_full)

$$\log f(x) = \log(\cos(x)) + x, \quad x \in ]-\pi/2, \pi/2[$, $r = $0.5, $\tau = $0.666, $r = $0.5, $\tau = $0.666, $r = $0.8, $\tau = $0.607, $r = $0.8, $\tau = $0.607, $r = $1.2, $\tau = $0.525, $r = $1.2, $\tau = $0.525, $r = $2, $\tau = $0.411, $r = $2, $\tau = $0.411$

Study of the acceptance ratio.¶

We study the acceptance ratio $\tau(r)$ as a function of the parameter $r$ for some distributions. In these cases, the value is exact and computed from the volume of $A_{f,r}$ and the one of $\tilde{A}_{f,r}$ . We draw the function $r \rightarrow \tau(r)$ . We conclude that the default value $r=1$ seems to be correct in all these cases.

mixt_gaus = ot.Mixture([ot.Normal(-2.0, 1.0), ot.Normal(4.0, 1.0)], [0.2, 0.8])
logPDF_list = [
    ot.Normal().getLogPDF(),
    ot.Normal(2).getLogPDF(),
    ot.Normal(3).getLogPDF(),
    ot.Exponential().getLogPDF(),
    mixt_gaus.getLogPDF(),
]
range_list = [
    ot.Normal().getRange(),
    ot.Normal(2).getRange(),
    ot.Normal(3).getRange(),
    ot.Exponential().getRange(),
    mixt_gaus.getRange(),
]
name_list = ["Normal(1)", "Normal(2)", "Normal(3)", "Exponential()", "Mixture Normal"]

g = ot.Graph(
    r"Acceptance ratio $\tau$ as a function of $r$", r"$r$", r"$\tau$", True, "topright"
)
g.setLogScale(1)
numberR = 50
a = (1e2) ** (2.0 / numberR)
data_r_ratio = ot.Sample(numberR + 1, 2)
isScaled = True
for k in range(len(logPDF_list)):
    ratio_Algo = otexp.RatioOfUniforms(logPDF_list[k], range_list[k], isScaled)
    for i in range(numberR + 1):
        r = a ** (i - 0.5 * numberR)
        ratio_Algo.setR(r)
        accept_Ratio = ratio_Algo.getAcceptanceRatio()
        data_r_ratio[i] = [r, accept_Ratio]
    draw = ot.Curve(data_r_ratio)
    draw.setLegend(name_list[k])
    g.add(draw)

isScaled = False
ratio_Algo = otexp.RatioOfUniforms(log_UnscaledPDF, range_PDF, isScaled)
for i in range(numberR + 1):
    r = a ** (i - 0.5 * numberR)
    ratio_Algo.setR(r)
    accept_Ratio = ratio_Algo.getAcceptanceRatio()
    data_r_ratio[i] = [r, accept_Ratio]
draw = ot.Curve(data_r_ratio)
draw.setLegend(
    r"$f(x) = \cos(x)e^x\mathbf{1}_{\left[-\frac{\pi}{2}, \frac{\pi}{2}\right]}(x)$"
)
g.add(draw)

g.setLegendPosition("upper left")
g.setLegendCorner([1.0, 1.0])
view = otv.View(g)

$Acceptance ratio $\tau$ as a function of $r$$

Exploration of the optimization problems.¶

We illustrate how to explore the results of the optimization problems that have been solved to define $\tilde{A}_{f,r}$ . We consider the example defined by the function:

$f(x) = \frac{1}{2} (2 + \sin(x)^2) \exp \left[- \left(2 + \cos(3x)^3 + \sin(2x)^3 \right) x \right] \mathbf{1}_{[0, 2 \pi]}(x).$

The method initialize returns all the starting points that have been used to compute respectively $u_{sup}$ , $v_{inf}$ and $v_{sup}$ . The maximum number of starting points is defined in ResourceMap, entry RatioOfUniformas-MaximumMultiStart. These points are selected scrolling through the Sobol low discrepance sequence in dimension $\inputDim = 1$ and which maximum size is defined in ResourceMap, entry RatioOfUniforms-CandidateNumber.

f = ot.SymbolicFunction(
    "x", "0.5 * (2 + sin(x)^2) * exp( -( 2 + cos(3*x)^3 + sin(2*x)^3) * x )"
)
log_UnscaledPDF = ot.ComposedFunction(ot.SymbolicFunction("x", "log(x)"), f)
lower_bound = 0.0
upper_bound = 2.0 * pi
range_PDF = ot.Interval(lower_bound, upper_bound)

graph = f.draw(lower_bound, upper_bound, 512)
graph.setTitle(
    r"Christian Robert function: $f(x) =  0.5(2 + \sin^2 x) e^{ -x( 2 + \cos^3 3x + \sin^3 2x)}, x \in [0, 2\pi]$"
)
graph.setXTitle(r"$x$")
graph.setYTitle(r"$f(x)$")
view = otv.View(graph)

ratio_Algo = otexp.RatioOfUniforms(log_UnscaledPDF, range_PDF)

$Christian Robert function: $f(x) = 0.5(2 + \sin^2 x) e^{ -x( 2 + \cos^3 3x + \sin^3 2x)}, x \in [0, 2\pi]$$

In order to get all the intermediate results, we change the ResourceMap, entry MultiStart-KeepResults.

Furthermore, as $I = [0, 2 \pi]= [a_i, b_i]$ , then $a_i \geq 0$ implies that $v_{inf} = 0$ . There are only the bounds $u_{sup}$ and $v_{sup}$ to find: the collection of optimization problems only contains 2 items.

ot.ResourceMap.SetAsBool("MultiStart-KeepResults", True)
coll_Multi_Start = ratio_Algo.initialize()
start_Sample_SupU = coll_Multi_Start[0].getStartingSample()
start_Sample_SupV = coll_Multi_Start[1].getStartingSample()

We plot on the previous graph all the starting points (red) used to find $u_{sup}$ and the associated maxima found by the algorithm. We also link each starting point to its associated optimal value (green).

Cloud_start_Usup = ot.Cloud(
    start_Sample_SupU, ot.Sample(start_Sample_SupU.getSize(), [0.0])
)
Cloud_start_Usup.setColor("red")
Cloud_start_Usup.setPointStyle("fsquare")
graph.add(Cloud_start_Usup)

coll_results = coll_Multi_Start[0].getResultCollection()
sample_optima = ot.Sample(0, 2)
for i in range(coll_results.getSize()):
    opt_x = coll_results[i].getOptimalPoint()[0]
    opt_value = f([opt_x])[0]
    sample_optima.add([opt_x, opt_value])
    line_convergence = ot.Curve([[start_Sample_SupU[i, 0], 0.0], [opt_x, opt_value]])
    line_convergence.setLineStyle("dashed")
    line_convergence.setColor("black")
    graph.add(line_convergence)
cloud_optima = ot.Cloud(sample_optima)
cloud_optima.setColor("green")
cloud_optima.setPointStyle("star")
graph.add(cloud_optima)
graph.setYTitle("")
view = otv.View(graph)

$Christian Robert function: $f(x) = 0.5(2 + \sin^2 x) e^{ -x( 2 + \cos^3 3x + \sin^3 2x)}, x \in [0, 2\pi]$$

view.ShowAll()

OpenTURNS

An Open source initiative for the Treatment of Uncertainties, Risks'N Statistics

Table of Contents

Previous topic

Next topic

This Page

Use the Ratio of Uniforms algorithm to sample a distribution¶

Study of the $(u,v)$ -space.¶

Study of the $(x,y)$ -space.¶

Illustration.¶

Study of the acceptance ratio.¶

Exploration of the optimization problems.¶

OpenTURNS

An Open source initiative for the Treatment of Uncertainties, Risks'N Statistics

Use the Ratio of Uniforms algorithm to sample a distribution¶

Study of the -space.¶

Study of the -space.¶

Illustration.¶

Study of the acceptance ratio.¶

Exploration of the optimization problems.¶

Study of the $(u,v)$ -space.¶

Study of the $(x,y)$ -space.¶