User-defined constants

Modeling
1

I have a working Stan program that includes Gauss-Legendre quadrature (GLQ) in a set of user-defined functions. GLQ requires a set of nodes and weights that never vary. I do not wish to store those values separately from my Stan program and provide them as data because I will be using it in all sorts of applications. And I certainly do not wish to recalculate them every time I run the program.

So far, using GLQ of order 20 (which needs 20 nodes and 20 weights) I laboriously typed them into the transformed data block. But I need to work with order of at least 64 and maybe much higher.

My question is whether there is an easy way to provide these constants to my Stan program.

If it helps, the program will be in an R package.

2

The best answer is that your package should probably be using the Gaussian-Legendre quadrature functions from Boost

https://www.boost.org/doc/libs/1_66_0/libs/math/doc/html/math_toolkit/gauss.html

via the external C++ mechanism

https://cran.r-project.org/package=rstan/vignettes/external.html

But if you just want something quickly, I would copy the values from a high-quality source, such as

http://www.holoborodko.com/pavel/numerical-methods/numerical-integration/#gauss_quadrature_abscissas_table

and then paste them into a spreadsheet. Delete the ±, add the commas, and paste back into a Stan program. Here it is for 64:

transformed data {
  real abscissae[32] = {
    0.024350292663425,
    0.072993121787799,
    0.121462819296121,
    0.169644420423993,
    0.217423643740007,
    0.264687162208767,
    0.311322871990211,
    0.357220158337668,
    0.402270157963992,
    0.446366017253464,
    0.489403145707053,
    0.531279464019895,
    0.571895646202634,
    0.611155355172393,
    0.648965471254657,
    0.685236313054233,
    0.719881850171611,
    0.752819907260532,
    0.783972358943341,
    0.813265315122798,
    0.84062929625258,
    0.865999398154093,
    0.889315445995114,
    0.910522137078503,
    0.92956917213194,
    0.946411374858403,
    0.961008799652054,
    0.973326827789911,
    0.983336253884626,
    0.991013371476745,
    0.996340116771955,
    0.999305041735772
  };
  real weights[32] = {
    0.04869095700914,
    0.048575467441504,
    0.048344762234803,
    0.047999388596458,
    0.04754016571483,
    0.04696818281621,
    0.046284796581315,
    0.045491627927418,
    0.044590558163757,
    0.043583724529324,
    0.042473515123654,
    0.041262563242624,
    0.03995374113272,
    0.038550153178616,
    0.03705512854024,
    0.035472213256882,
    0.033805161837142,
    0.032057928354852,
    0.030234657072403,
    0.02833967261426,
    0.026377469715055,
    0.024352702568711,
    0.022270173808383,
    0.02013482315353,
    0.017951715775697,
    0.015726030476025,
    0.013463047896719,
    0.011168139460131,
    0.008846759826364,
    0.006504457968978,
    0.004147033260562,
    0.001783280721696
  };
}
1 Like
3

Excellent answer. Many thanks, Ben.