Skip to content

Commit

Permalink
Reverted to truncated normal example and cleaned up latex (Issue #2466)
Browse files Browse the repository at this point in the history
  • Loading branch information
@bbbales2
bbbales2 committed Sep 14, 2018
1 parent 52c2590 commit 4a0595a
Show file tree
Hide file tree
Showing 2 changed files with 79 additions and 78 deletions.
14 changes: 0 additions & 14 deletions src/docs/bibtex/all.bib
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1356,17 +1356,3 @@ @article{Tanaka:2009
doi="10.1007/s00211-008-0195-1",
url="https://doi.org/10.1007/s00211-008-0195-1"
}

@article{Rothwell:2008,
author = {Rothwell, E. J.},
title = {Computation of the logarithm of Bessel functions of complex argument and fractional order},
journal = {Communications in Numerical Methods in Engineering},
volume = {24},
number = {3},
pages = {237-249},
keywords = {Bessel functions, numerical analysis, electromagnetic scattering, boundary value problems, eigenvalues and eigenfunctions},
doi = {10.1002/cnm.972},
url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/cnm.972},
eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1002/cnm.972},
abstract = {Abstract Algorithms for computing the logarithm of Bessel functions of complex argument and fractional order are presented. These algorithms will be useful when direct evaluation of the Bessel function causes computer overflow or underflow. Several numerical results are presented, and an example of electromagnetic scattering by a conducting wedge is given. Copyright © 2006 John Wiley \& Sons, Ltd.}
}
143 changes: 79 additions & 64 deletions src/docs/users-guide/examples.tex
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8934,47 +8934,28 @@ \subsection{Maximum Number of Steps}
\chapter{Computing One Dimensional Integrals}

\noindent
Definite and indefinite one dimensional integrals can be performed in Stan via the \code{integrate\_1d} function.
Definite and indefinite one dimensional integrals can be performed in Stan using the \code{integrate\_1d} function.

As an example, the Mat\'ern covariance kernel can be written
As an example, the normalizing constant of a left-truncated normal distribution is

\begin{align}
k_\nu(r) = \frac{2^{1 - \nu}}{\Gamma(\nu)} \left( \frac{\sqrt{2 \nu} r}{l} \right)^\nu K_\nu \left( \frac{\sqrt{2 \nu} r}{l} \right)
\end{align}

\noindent where K_\nu is a modified Bessel function of the second time. For general cases of $\nu$, this Bessel function can be computed
with the definite integral given in \cite{Rothwell:2008}.

\begin{align}
K_\nu(z) = \frac{\pi}{\Gamma(\nu + \frac{1}{2})} (2 z)^{-\nu} e^{-z} \int_0^1 [\beta e^{-x^\beta} (2 z + x^\beta)^{\nu - \frac{1}{2}} x^{n - 1} e^{-\frac{1}{x}} x^{-2 \nu - 1} (2 z x + 1)^{\nu - \frac{1}{2}}] dx
\end{align}
\begin{align*}
\int_a^\infty \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-\frac{1}{2}\frac{(x - \mu)^2}{\sigma^2}}
\end{align*}

\noindent where
\noindent
To compute this integral in Stan, the integrand must first be defined as a Stan function (see \refchapter{functions-programming} for more information on coding user-defined functions).

\begin{align}
\beta = \frac{2 n}{2 \nu + 1}
\end{align}
\begin{stancode}
real normal_density(real x, // Function argument
real xc, // Complement of function argument
// on the domain (defined later)
real[] theta, // parameters
real[] x_r, // data (real)
int[] x_i) { // data (integer)
real mu = theta[1];
real sigma = theta[2];

\noindent
$n$ is set to $8$, mimicking \ref{Rothwell:2008}. To compute this integral in Stan, the integrand must first be defined as a Stan function (see \refchapter{functions-programming} for more information on coding user-defined functions).

\begin{stancode}
real Kv_int(real x, // Function argument
real xc, // Complement of function argument
// on the domain (defined later)
real[] theta, // parameters
real[] x_r, // data (real)
int[] x_i) { // data (integer)
real nu = theta[1];
real beta = 16.0 / (2.0 * nu + 1.0);

real first = beta * exp(-pow(x, beta)) * pow(2 * z + pow(x, beta), nu - 0.5) * pow(x, 7);
real second = exp(-1.0 / x);
if (second > 0) {
second = second * pow(x, -2 * nu - 1) * pow(2 * z * x + 1, v - 0.5);
}

return first + second;
return 1 / (sqrt(2 * pi()) * sigma) * exp(-0.5 * ((x - mu) / sigma)^2);
}
\end{stancode}
%
Expand All @@ -8986,7 +8967,8 @@ \chapter{Computing One Dimensional Integrals}
(see \refsection{integral-precision} for details on how to use this). The argument
\code{theta} is used to pass in arguments of the integral
that are a function of the parameters in our model. The arguments \code{x\_r}
and \code{x\_i} are likewise used to pass in arguments of the integral that are data.
and \code{x\_i} are used to pass in real and integer arguments of the integral that are
not a function of our parameters.

\subsubsection{Strict Signature}
%
Expand All @@ -8997,7 +8979,7 @@ \subsubsection{Strict Signature}
for data or a zero-length vector for parameters if the integral does not involve
data or parameters.

\section{Calling the 1D Integrator}
\section{Calling the Integrator}
%
Suppose that our model requires evaluating the lpdf of a left-truncated normal, but
the truncation limit is to be estimated as a parameter. Because the truncation
Expand All @@ -9006,7 +8988,8 @@ \section{Calling the 1D Integrator}
1D integrator. The more efficient way to perform the correct normalization in Stan
is described in \refchapter{truncated-data}.

Such a model might look like (including the function defined above):
Such a model might look like (include the function defined at the beginning of this
chapter to make this code compile):

%
\begin{stancode}
Expand All @@ -9031,18 +9014,22 @@ \section{Calling the 1D Integrator}
sigma ~ normal(0, 1);
left_limit ~ normal(0, 1);
target += normal_lpdf(y | mu, sigma);
target += log(integrate_1d(normal_density, left_limit, positive_infinity(), { mu, sigma }, x_r, x_i));
target += log(integrate_1d(normal_density,
left_limit,
positive_infinity(),
{ mu, sigma }, x_r, x_i));
}
\end{stancode}
%
\subsection{Limits of integration}
The limits of integration can be finite or infinite. The infinite limits are made available via
the Stan calls \code{negative\_infinity()} and \code{positive\_infinity()}.
\subsection{Limits of Integration}
The limits of integration can be finite or infinite. The infinite limits are
made available via the Stan calls \code{negative\_infinity()} and
\code{positive\_infinity()}.

If both limits are either \code{negative\_infinity()} or \code{positive\_infinity()}, the integral
and its gradients are set to zero.
If both limits are either \code{negative\_infinity()} or
\code{positive\_infinity()}, the integral and its gradients are set to zero.

\subsection{Data versus Parameters}
\subsection{Data Versus Parameters}
%
The arguments for the real data \code{x\_r} and the integer data \code{x\_i}
must be expressions that only involve data or transformed data variables.
Expand All @@ -9053,35 +9040,50 @@ \subsection{Data versus Parameters}
derivatives of the integral with respect to the endpoints are handled
with the Leibniz integral rule).

\section{Integrator convergence}
\section{Integrator Convergence}
%
The integral is performed with the iterative 1D quadrature methods implemented in the Boost library \cite{BoostQuadrature:2017}. If the nth estimate of the integral is denoted $I_n$ and the nth estimate of the norm of the integral is denoted $|I|_n$, the iteration is terminated when
The integral is performed with the iterative 1D quadrature methods implemented
in the Boost library \cite{BoostQuadrature:2017}. If the nth estimate of the
integral is denoted $I_n$ and the nth estimate of the norm of the integral is
denoted $|I|_n$, the iteration is terminated when

\begin{align*}
\frac{{|I_{n + 1} - I_n|}}{{|I|_{n + 1}}} < \text{relative tolerance}
\end{align*}

\noindent
The \code{relative\_tolerance} parameter can be optionally specified as the last argument to \code{integrate\_1d} and defaults to the square root of the machine epsilon of double precision floating point numbers (that is a relative tolerance of around \code{1.5e-8}).
The \code{relative\_tolerance} parameter can be optionally specified as the
last argument to \code{integrate\_1d} and defaults to the square root of the
machine epsilon of double precision floating point numbers (that is a relative
tolerance of around \code{1.5e-8}).

\subsection{Zero-crossing integrals} \label{zero-crossing.section}
Integrals on the (possibly infinite) interval $(a, b)$ that cross zero are split into two integrals, one from $(a, 0)$ and one from $(0, b)$. This is because the quadrature methods employed internally can have difficulty near zero.
\subsection{Zero-crossing Integrals} \label{zero-crossing.section}
Integrals on the (possibly infinite) interval $(a, b)$ that cross zero are
split into two integrals, one from $(a, 0)$ and one from $(0, b)$. This is
because the quadrature methods employed internally can have difficulty near
zero.

In this case, each integral is separately integrated to the given \code{relative\_tolerance}.
In this case, each integral is separately integrated to the given
\code{relative\_tolerance}.

\subsection{Avoiding precision loss near limits of integration in definite integrals} \label{integral-precision.section}
\subsection{Avoiding precision loss near limits of integration in definite
integrals} \label{integral-precision.section}

If care is not taken, the quadrature can suffer from numerical loss of precision near the endpoints of definite integrals.
If care is not taken, the quadrature can suffer from numerical loss of
precision near the endpoints of definite integrals.

For instance, in integrating the pdf of a beta distribution when the values of $\alpha$ and $\beta$ are small, most of the probability mass is lumped near zero and one.
For instance, in integrating the pdf of a beta distribution when the values of
$\alpha$ and $\beta$ are small, most of the probability mass is lumped near zero
and one.

The pdf of a beta distribution is proportional to
\begin{align*}
p(x) \propto x^{\alpha - 1}(1 - x)^{\beta - 1}
\end{align*}

\noindent
Normalizing this distribution requires computing the integral of $p(x)$ from zero to one. In Stan code, the integrand might look like:
Normalizing this distribution requires computing the integral of $p(x)$ from
zero to one. In Stan code, the integrand might look like:

%
\begin{stancode}
Expand All @@ -9095,9 +9097,15 @@ \subsection{Avoiding precision loss near limits of integration in definite integ
%

\noindent
The issue is that there will be numerical breakdown in the precision of \code{1.0 - x} as \code{x} gets close to one. This is because of the limited precision of double precision floating numbers. This integral will fail to converge for values of \code{alpha} and \code{beta} much less than one.
The issue is that there will be numerical breakdown in the precision of
\code{1.0 - x} as \code{x} gets close to one. This is because of the limited
precision of double precision floating numbers. This integral will fail to
converge for values of \code{alpha} and \code{beta} much less than one.

This is where \code{xc} is useful. It is defined, for definite integrals, as a high precision version of the distance from \code{x} to the nearest endpoint. To make use of this for the beta integral, the integrand can be re-coded:
This is where \code{xc} is useful. It is defined, for definite integrals,
as a high precision version of the distance from \code{x} to the nearest
endpoint. To make use of this for the beta integral, the integrand can be
re-coded:

%
\begin{stancode}
Expand All @@ -9118,9 +9126,16 @@ \subsection{Avoiding precision loss near limits of integration in definite integ
%

\noindent
This version of the integrand will converge for much smaller values of \code{alpha} and \code{beta} than otherwise possible.

Note, \code{xc} is only used for definite integrals. If either the left endpoint is at negative infinity or the right endpoint is at positive infinity, \code{xc} will be NaN.

For zero-crossing definite integrals (see \refsection{zero-crossing}) the integrals are broken into two pieces ($(a, 0)$ and $(0, b)$ for endpoints $a < 0$ and $b > 0$) and \code{xc} is a high precision version of the distance to the limits of each of the two integrals separately. This means \code{xc} will be the a high precision version of \code{a - x}, \code{x}, or \code{b - x}, depending on the value of x and the endpoints.

This version of the integrand will converge for much smaller values of
\code{alpha} and \code{beta} than otherwise possible.

Note, \code{xc} is only used for definite integrals. If either the left endpoint
is at negative infinity or the right endpoint is at positive infinity, \code{xc}
will be NaN.

For zero-crossing definite integrals (see \refsection{zero-crossing}) the
integrals are broken into two pieces ($(a, 0)$ and $(0, b)$ for endpoints
$a < 0$ and $b > 0$) and \code{xc} is a high precision version of the distance
to the limits of each of the two integrals separately. This means \code{xc} will
be the a high precision version of \code{a - x}, \code{x}, or \code{b - x},
depending on the value of x and the endpoints.

0 comments on commit 4a0595a

Please sign in to comment.