lintegrate

一个数值积分库，用于您需要在对数空间中处理需要积分的函数的自然对数时。

此库提供了三个数值积分函数，基于GSL函数，当只有函数的自然对数给出时，可以对其进行积分，并返回该积分的自然对数。这三个函数

• lintegrate_qag
• lintegrate_qng
• lintegrate_cquad

分别是GSL函数

• gsl_integration_qag
• gsl_integration_qng
• gsl_integration_cquad

的等价函数。这些函数在例如，您可以计算高斯似然函数的自然对数（在高斯函数的指数会导致零或无穷大时）时非常有用，并且您希望数值地找到高斯函数本身的积分。

lintegrate_qag、lintegrate_qng和lintegrate_cquad函数都有包装函数（在名称后附加_split），允许用户指定在计算时将积分分割成的一组区间。这些区间可以是对数空间中均匀分布的，对于积分函数在接近零时具有非常明显的峰值的情况非常有用。

完整的API文档和示例可以在此处找到。

示例

使用这些函数的示例如下

/* example using lintegrate functionality */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>

#include <lintegrate.h>

/* create function for integration */
double lintegrand(double x, void *params);

struct intparams {
  double mu;
  double sig;
};

double lintegrand(double x, void *params){
  struct intparams * p = (struct intparams *)params;
  double mu = p->mu;
  double sig = p->sig;

  return -0.5*(mu-x)*(mu-x)/(sig*sig);
}

double integrand(double x, void *params){
  struct intparams * p = (struct intparams *)params;
  double mu = p->mu;
  double sig = p->sig;

  return exp(-0.5*(mu-x)*(mu-x)/(sig*sig));
}

int main( int argv, char **argc ){
  gsl_function F;
  struct intparams params;
  gsl_integration_workspace *w = gsl_integration_workspace_alloc (100);
  gsl_integration_cquad_workspace *cw = gsl_integration_cquad_workspace_alloc(50);
  double qaganswer = 0., qnganswer = 0., cquadanswer = 0., answer = 0.;
  double abserr = 0.;
  size_t neval = 0;

  double minlim = -6.; /* minimum for integration range */
  double maxlim = 6.;  /* maximum for integration range */

  double abstol = 1e-10; /* absolute tolerance */
  double reltol = 1e-10; /* relative tolerance */

  params.mu = 0.;
  params.sig = 1.;

  F.function = &lintegrand;
  F.params = &params;

  /* integrate log of function using QAG */
  lintegration_qag(&F, minlim, maxlim, abstol, reltol, 100, GSL_INTEG_GAUSS31, w, &qaganswer, &abserr);

  /* integrate log of function using QNG */
  lintegration_qng(&F, minlim, maxlim, abstol, reltol, &qnganswer, &abserr, &neval);

  /* integrate log of function using CQUAD */
  lintegration_cquad(&F, minlim, maxlim, abstol, reltol, cw, &cquadanswer, &abserr, &neval);

  /* integrate function using GSL QAG */
  F.function = &integrand;
  gsl_integration_qag(&F, minlim, maxlim, abstol, reltol, 100, GSL_INTEG_GAUSS31, w, &answer, &abserr);

  gsl_integration_workspace_free(w);
  gsl_integration_cquad_workspace_free(cw);

  fprintf(stdout, "Answer \"lintegrate QAG\" = %.8lf\n", qaganswer);
  fprintf(stdout, "Answer \"lintegrate QNG\" = %.8lf\n", qnganswer);
  fprintf(stdout, "Answer \"lintegrate CQUAD\" = %.8lf\n", cquadanswer);
  fprintf(stdout, "Answer \"gsl_integrate_qag\" = %.8lf\n", log(answer));
  fprintf(stdout, "Analytical answer = %.8lf\n", log(sqrt(2.*M_PI)));

  return 0;
}

要求

• GSL - 在Debian/Ubuntu (16.04)上使用以下命令安装： sudo apt-get install libgsl-dev

安装

可以通过在基本目录中输入 sudo scons 使用 scons 构建库。要全局安装库（默认情况下在 /usr/local/lib），请运行

sudo scons
sudo scons install

用户可以从源代码构建和安装包含函数包装器的Python模块，例如

pip install .

在存储库目录中。

Python模块也可以通过pip从PyPI安装

pip install lintegrate

或在Conda环境中使用

conda install -c conda-forge lintegrate

Python

如果已安装Python模块，则具有以下函数

• lqng - 对 lintegration_qng 的包装
• lqag - 对 lintegration_qag 的包装
• lcquad - 对 lintegration_cquad 的包装
• logtrapz - 在值网格上使用梯形规则进行积分

lqng、lqag 和 lcquad 函数的使用方式与scipy的 quad 函数类似。

它们的使用示例可以是

from lintegrate import lqag, lqng, lcquad, logtrapz
import numpy as np

# define the log of the function to be integrated
def integrand(x, args):
    mu, sig = args # unpack extra arguments
    return -0.5*((x-mu)/sig)**2

# set integration limits
xmin = -6.
xmax = 6.

# set additional arguments
mu = 0.
sig = 1.

resqag = lqag(integrand, xmin, xmax, args=(mu, sig))
resqng = lqng(integrand, xmin, xmax, args=(mu, sig))
rescquad = lcquad(integrand, xmin, xmax, args=(mu, sig))
restrapz = logtrapz(integrand, np.linspace(xmin, xmax, 100), args=(mu, sig))

R

在 R 中，可以使用 reticulate 包调用 lintegrate 中的函数。上面的示例将是

library(reticulate)
py_install("lintegrate", pip = TRUE) ## run once to make sure lintegrate is installed and visible to reticulate.
lint <- import("lintegrate", convert = FALSE)
integrand <- function(x, args){
  mu = args[1]
  sig = args[2]
  return(-.5 * ((x-mu)/sig)^2 )
} 
integrand <- Vectorize(integrand)
mu <- 0
sig <- 1
mmin <- -10
mmax <- 10
lint$lqag(py_func(integrand), r_to_py(mmin), r_to_py(mmax), c(mu, sig))

引用

如果在研究中使用 lintegrate，我将不胜感激，如果您引用了与软件相关的 JOSS 论文。以下BibTeX引用可以使用

@article{Pitkin2022,
  doi = {10.21105/joss.04231},
  url = {https://doi.org/10.21105/joss.04231},
  year = {2022},
  publisher = {The Open Journal},
  volume = {7},
  number = {73},
  pages = {4231},
  author = {Matthew Pitkin},
  title = {lintegrate: A C/Python numerical integration library for working in log-space},
  journal = {Journal of Open Source Software}
}

您还可以引用GSL参考《M. Galassi等人，GNU科学库参考手册（第3版），ISBN 0954612078》和URL https://gnu.ac.cn/software/gsl/。

© 2017 Matthew Pitkin