时间戳

首页 > 解决方案 > `integrate_1d` 语法并以交互方式调试它

integration - `integrate_1d` 语法并以交互方式调试它

问题描述

我是第一次尝试设置自定义可能性并使用integrate_1d. 根据文档，我将我的被积函数设置如下：

// internal function to integrate over z 
// integrate_1d is super picky: the argument types need to be real, real, 
    real[], real[], int[]
// e.g., cannot use vectors in place of real[] or int[]

real integrand( real x,  // the observed Z-stat; value at which to evaluate integral
              real xc,
              real[] theta,
              real[] x_r,
              int[] x_i ){

    // separate the parameters
    real zeta = theta[1];
    real eta = theta[2];
    real vi = theta[3];

    // significance indicator
    int signif; 
    signif = fabs(x) > 1.96;

    return( (1/eta) * (1 - signif) + signif ) * exp( normal_lpdf( x | zeta, sqrt(vi) ) );
  }

我正在尝试测试和调试此功能，因为它最终会导致在我运行模型时评估对数似然性时出现错误。为此，我希望能够以交互方式运行该函数，传递我自己选择的参数。我正在尝试使用 rstanexpose_stan_functions来做到这一点：

code = "functions{

  real integrand( real x,  // the observed Z-stat
                  real xc,
                  real[] theta,
                  real[] x_r,
                  int[] x_i ){

    // separate the parameters
    real zeta = theta[1];
    real eta = theta[2];
    real vi = theta[3];

    // significance indicator
    int signif; 
    signif = fabs(x) > 1.96;

    return( (1/eta) * (1 - signif) + signif ) * exp( normal_lpdf( x | zeta, sqrt(vi) ) );
  }
}"

expose_stan_functions(stanc(model_code = code))

integrand(1, {1,1,1}, real[], int[], 1e-8)

integrate_1d( integrand, negative_infinity(), positive_infinity(), {1, 1, 1}, real[], int[], 1e-8 )

但是调用integrand和integrate_1d每次都抛出错误

错误：意外的“，”

就好像我用根本不正确的语法调用函数一样。（在调用 时integrate_1d，我不太确定如何处理xc。我省略了它，因为文档（第 9.3.2.4 节）似乎表明它应该是NaN一个不定积分。）

integrand如果expose_stan_functions不是理想的方式，我也欢迎完全不同的交互式调试建议。

标签： integrationstanrstan

解决方案

该integrate_1d函数是一个 Stan / C++ 函数，所以如果你想从 R 调用它，你需要指定另一个调用它的函数。在你的情况下，你可以做

code = "functions{

  real integrand( real x,  // the observed Z-stat
                  real xc,
                  real[] theta,
                  real[] x_r,
                  int[] x_i ){

    // separate the parameters
    real zeta = theta[1];
    real eta = theta[2];
    real vi = theta[3];

    // significance indicator
    int signif; 
    signif = fabs(x) > 1.96;

    return( (1/eta) * (1 - signif) + signif ) * exp( normal_lpdf( x | zeta, sqrt(vi) ) );
  }

  real area(real[] theta, data real[] x_r) {
    int x_i[0];
    return integrate_1d(integrand, negative_infinity(), positive_infinity(),
                        theta, x_r, x_i, 1e-8);
  }
}"

library(rstan)
expose_stan_functions(stanc(model_code = code))
area(theta = rep(1, 3), x_r = double()) # yields 1

推荐阅读