What is Gaussian quadrature and how to implement it in Stan?
Happy New "Reiwa" Era!
Numerical integral is necessary to do analysis with survival models such as joint modeling, however it is difficult to find relevant information. Then, I write some tips on that. Don't be blind on that!
Suppose the function \(f(x)\) can be defined as the following 3th order polynomial function: \begin{align} f(x):= a_0 + a_1 x + a_2 x^2 + a_3 x^3, \nonumber \end{align} and its analytical solution for the integral with the interval \([-1,1]\) is: \begin{align} \int_{-1}^1 f(x) dx &= \int_{-1}^1 a_0 + a_1 x + a_2 x^2 + a_3 x^3 dx \nonumber \\ &= \left[ a_0 x + \frac{1}{2}a_1 x^2 + \frac{1}{3}a_2 x^3 + \frac{1}{4}a_3 x^4 \right]_{-1}^1 \nonumber \\ &= 2 a_0 + \frac{2}{3}a_2 \label{eq:approx} \end{align} On the other hand, Eq.(\ref{eq:c1c2}) yields \begin{align} \int_{-1}^1 f(x) dx &\approx C_1 \Bigl( a_0 + a_1 x_1 + a_2 x_1^2 + a_3 x_1^3 \Bigr) + C_2 \Bigl( a_0 + a_1 x_2 + a_2 x_2^2 + a_3 x_2^3 \Bigr) \nonumber \\ &= \Bigl( C_1 + C_2 \Bigr) a_0 + \Bigl( C_1 x_1 + C_2 x_2 \Bigr) a_1 \nonumber \\ &\ \ \ \ \ \ + \Bigl( C_1 x_1^2 + C_2 x_2^2 \Bigr) a_2 + \Bigl( C_1 x_1^3 + C_2 x_2^3 \Bigr) a_3. \label{eq:c1c1c2c2} \end{align} By comparing Eq.(\ref{eq:approx}) and Eq.(\ref{eq:c1c1c2c2}), we get the following 4 equalities. \begin{align} 2 &=C_1 + C_2 \nonumber \\ 0 &= C_1 x_1 + C_2 x_2 \nonumber \\ \frac{2}{3} &= C_1 x_1^2 + C_2 x_2^2 \nonumber \\ 0 &= C_1 x_1^3 + C_2 x_2^3 \nonumber \end{align} Although it might have multiple solutions, assume \(C_1=C_2\) and \(x_1 \le x_2\), \begin{align} 2 =C_1 + C_2 \ \ \ \Longleftrightarrow \ \ \ 2 =2C_1 \ \ \ \Longleftrightarrow \ \ \ C_1=C_2=1,\nonumber \\ \end{align} and then \begin{align} 0 = C_1 x_1 + C_2 x_2 \ \ \ \Longleftrightarrow \ \ \ 0 = x_1 + x_2 \ \ \ \Longleftrightarrow \ \ \ x_2=-x_1,\nonumber \\ \end{align} and \begin{align} \frac{2}{3} = C_1 x_1^2 + C_2 x_2^2 \ \ \ \Longleftrightarrow \ \ \ \frac{2}{3} = 2x_1^2 \ \ \ \Longleftrightarrow \ \ \ \left\{ \begin{array}{l} x_1=-\sqrt{\frac{1}{3}} \\ x_2=+\sqrt{\frac{1}{3}} \end{array} \right.. \nonumber \\ \end{align} Then the integral is approximated by \begin{align} \int_{-1}^1 f(x) dx \approx f\left( -\sqrt{\frac{1}{3}} \right) + f\left( \sqrt{\frac{1}{3}} \right) \end{align} This is the 2-points Gaussian quadrature. in short, one can consider n-points Gaussian quadrature by assuming the function can be defined as \(2n-1\) order polynomial function. As in the case of 2-points Gaussian quadrature above, one can specify the weights and argument points with similar simple assumptions such as \(C_1=C_2\) and \(x_1 \ge x_2\) with 2-points case. The solutions up to 5-points are provided on the page at Wikipedia.
Numerical integral is necessary to do analysis with survival models such as joint modeling, however it is difficult to find relevant information. Then, I write some tips on that. Don't be blind on that!
Basic ideas for Gaussian quadrature
Consider the following integral: \begin{align} \int_{-1}^1 f(x) dx,\nonumber \end{align} where \(f\) is any function of \(x\). And let the function be able to be approximated by \begin{align} \int_{-1}^1 f(x) dx \approx C_1 f(x_1) + C_2 f(x_2) \label{eq:c1c2} \end{align} where \(C_1, C_2\) are some wights and \(x_1 (-1 \le x_1 \le 1),x_2 (-1 \le x_2 \le 1)\) are arguments which are representative points for the integral. And it also should be noted that we can change the upper and lower bound of the integral \([a,b]\) to \([-1,1]\) via the basic variable exchange rule. Suppose \(y(x)=\displaystyle \frac{b-a}{2}x + \frac{b+a}{2}\), \begin{align} \int_a^b f(x) dx &= \int_{y(a)}^{y(b)} f(y) \left| \frac{dy}{dx} \right| dx \nonumber \\ &= \left| \frac{d}{dx}\left( \frac{b-a}{2}x + \frac{b+a}{2} \right) \right|\int_{-1}^1 f \left( \frac{b-a}{2}x + \frac{b+a}{2} \right) dx\nonumber \\ &=\frac{b-a}{2} \int_{-1}^1 f \left( \frac{b-a}{2}x + \frac{b+a}{2} \right) dx \label{eq:valex} \end{align} Then what we need to consider is only the integral with the interval \([-1,1]\).Suppose the function \(f(x)\) can be defined as the following 3th order polynomial function: \begin{align} f(x):= a_0 + a_1 x + a_2 x^2 + a_3 x^3, \nonumber \end{align} and its analytical solution for the integral with the interval \([-1,1]\) is: \begin{align} \int_{-1}^1 f(x) dx &= \int_{-1}^1 a_0 + a_1 x + a_2 x^2 + a_3 x^3 dx \nonumber \\ &= \left[ a_0 x + \frac{1}{2}a_1 x^2 + \frac{1}{3}a_2 x^3 + \frac{1}{4}a_3 x^4 \right]_{-1}^1 \nonumber \\ &= 2 a_0 + \frac{2}{3}a_2 \label{eq:approx} \end{align} On the other hand, Eq.(\ref{eq:c1c2}) yields \begin{align} \int_{-1}^1 f(x) dx &\approx C_1 \Bigl( a_0 + a_1 x_1 + a_2 x_1^2 + a_3 x_1^3 \Bigr) + C_2 \Bigl( a_0 + a_1 x_2 + a_2 x_2^2 + a_3 x_2^3 \Bigr) \nonumber \\ &= \Bigl( C_1 + C_2 \Bigr) a_0 + \Bigl( C_1 x_1 + C_2 x_2 \Bigr) a_1 \nonumber \\ &\ \ \ \ \ \ + \Bigl( C_1 x_1^2 + C_2 x_2^2 \Bigr) a_2 + \Bigl( C_1 x_1^3 + C_2 x_2^3 \Bigr) a_3. \label{eq:c1c1c2c2} \end{align} By comparing Eq.(\ref{eq:approx}) and Eq.(\ref{eq:c1c1c2c2}), we get the following 4 equalities. \begin{align} 2 &=C_1 + C_2 \nonumber \\ 0 &= C_1 x_1 + C_2 x_2 \nonumber \\ \frac{2}{3} &= C_1 x_1^2 + C_2 x_2^2 \nonumber \\ 0 &= C_1 x_1^3 + C_2 x_2^3 \nonumber \end{align} Although it might have multiple solutions, assume \(C_1=C_2\) and \(x_1 \le x_2\), \begin{align} 2 =C_1 + C_2 \ \ \ \Longleftrightarrow \ \ \ 2 =2C_1 \ \ \ \Longleftrightarrow \ \ \ C_1=C_2=1,\nonumber \\ \end{align} and then \begin{align} 0 = C_1 x_1 + C_2 x_2 \ \ \ \Longleftrightarrow \ \ \ 0 = x_1 + x_2 \ \ \ \Longleftrightarrow \ \ \ x_2=-x_1,\nonumber \\ \end{align} and \begin{align} \frac{2}{3} = C_1 x_1^2 + C_2 x_2^2 \ \ \ \Longleftrightarrow \ \ \ \frac{2}{3} = 2x_1^2 \ \ \ \Longleftrightarrow \ \ \ \left\{ \begin{array}{l} x_1=-\sqrt{\frac{1}{3}} \\ x_2=+\sqrt{\frac{1}{3}} \end{array} \right.. \nonumber \\ \end{align} Then the integral is approximated by \begin{align} \int_{-1}^1 f(x) dx \approx f\left( -\sqrt{\frac{1}{3}} \right) + f\left( \sqrt{\frac{1}{3}} \right) \end{align} This is the 2-points Gaussian quadrature. in short, one can consider n-points Gaussian quadrature by assuming the function can be defined as \(2n-1\) order polynomial function. As in the case of 2-points Gaussian quadrature above, one can specify the weights and argument points with similar simple assumptions such as \(C_1=C_2\) and \(x_1 \ge x_2\) with 2-points case. The solutions up to 5-points are provided on the page at Wikipedia.
How to implement it in Stan
This is the easiest implementation in Stan, only core part was described in the following code. Suppose we have survival time \(t\) comprising \(n\) observations, and consider the case where an integral \(\int_0^t f(u)du\) should be included as log probability density, i.e., the survival probability with a hazard function \(\displaystyle \log\left(S(t)\right)=\log \left[ \exp \left(-\int_0^t h(u) du \right) \right] = -\int_0^t h(u) du\), and it is going to be approximated by 2-points Gaussian quadrature. (I know that 2-points Gaussian quadrature is not realistic, but it would be good as a toy example.) One can prepare the weights and arguments outside Rstan as follows:d = list( n = n, # the number of observations t = t, # the survival time C = c(1,1), # vector for weights X = c(-1/sqrt(3),-1/sqrt(3)) # vector for arguments )On the other hand, those variables can be used in Stan as shown below. The part \(\displaystyle \frac{b-a}{2}\) and \(\displaystyle \frac{b+a}{2}\) in Eq.(\ref{eq:valex}) become both \(\displaystyle \frac{t_i \pm 0}{2}=\frac{t_i}{2}\) as shown in the lines 14 and 15. Note that the function f in the line 20 must be either defined beforehand or written directly on that part.
model_code='
data {
int n;
real t[n];
real C[2];
real X[2];
}
transformed parameters{
real X1[n];
real X2[n];
for (i in 1:n)
X1[i] = t[i]/2*X[1]+t[i]/2; // variable exchange
X2[i] = t[i]/2*X[2]+t[i]/2; // variable exchange
}
model {
for (i in 1:n)
target += t[i]/2*(C[1]*f(X1[i])+C[2]*f(X2[i])); // Eq.(2)
}
'
It can be applied to n-points Gaussian quadrature and n-points Gauss-Kronrod Quadrature. The weights and arguments are provided on the site. Should keep it mind that it takes a long time to execute as the number of points increases.
5-points Gaussian quadrature
Inmost cases, especially for model development stage, 5-points Gaussian quadrature gives enough close approximation for the integration. You can use the following code if necessary.X = c( 1/3*sqrt(5 - 2*sqrt(10/7)), -1/3*sqrt(5 - 2*sqrt(10/7)), 0, 1/3*sqrt(5 + 2*sqrt(10/7)), -1/3*sqrt(5 + 2*sqrt(10/7)) ) C =c( ((322+13*sqrt(70))/900), ((322+13*sqrt(70))/900), (128/225), ((322-13*sqrt(70))/900), ((322-13*sqrt(70))/900) )
7-point Gauss rule with a 15-point Kronrod
If you need to use more precise approximation, you can use 7-point Gauss rule with a 15-point Kronrod, the numbers below are from WikipediaX = c( -0.991455371120813, -0.949107912342759, -0.864864423359769, -0.741531185599394, -0.586087235467691, -0.405845151377397, -0.207784955007898, 0, 0.207784955007898, 0.405845151377397, 0.586087235467691, 0.741531185599394, 0.864864423359769, 0.949107912342759, 0.991455371120813 ) C =c( 0.022935322010529, 0.063092092629979, 0.104790010322250, 0.140653259715525, 0.169004726639267, 0.190350578064785, 0.204432940075298, 0.209482141084728, 0.204432940075298, 0.190350578064785, 0.169004726639267, 0.140653259715525, 0.104790010322250, 0.063092092629979, 0.022935322010529 )
Comments
Post a Comment