I am trying to write a custom distribution function (lpdf) in C++ that is called from R with rstan. However, when I compile and fit it, I am getting some errors. Here the contents of my files.
Rscript:
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
model_code = "
functions {
real custom_function_lpdf(real y, real A, real B);
}
data {
int N;
real y[N];
}
parameters {
real A;
real B;
}
model {
for(i in 1:N){
y[i] ~ custom_function(A, B);
}
}
"
data <- list(N = 1, A = 1, B = 1)
model <- stan_model(model_code = model_code,
allow_undefined = TRUE,
includes = paste0('\n#include "', file.path(getwd(), 'custom_function_lpdf.hpp'), '"\n'))
fit <- sampling(model, data = data, iter=1000, chains=1)
print(fit)
custom_function_lpdf.hpp:
#include <boost/math/tools/promotion.hpp>
using namespace boost::math;
using namespace std;
template <typename T_y, typename T_A, typename T_B>
typename boost::math::tools::promote_args<T_y, T_A, T_B>::type
parabola_function_lpdf(const T_y& y,
const T_A& A,
const T_B& B,
std::ostream* pstream__) {
return A*B*y;
}
I get the following error:
Error in dyn.load(libLFile) :
unable to load shared object '/tmp/RtmpKiWiMx/file2dee16490e35f8.so':
/tmp/RtmpKiWiMx/file2dee16490e35f8.so: undefined symbol: _ZN62model2dee1627232579_76aac3ff274b7f875f40dd19472b96f5_namespace20custom_function_lpdfILb0EdddEEN5boost4math5tools12promote_argsIT0_T1_T2_fffE4typeERKS5_RKS6_RKS7_PSo
The C++ code seems to compile but the error is related to linking. What am I doing wrong? Is there any working example of external custom C++ distribution function?
For example, if I want to implement the following std_normal_lpdf function from Stan as a custom_function_lpdf, how can I do it?
#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/constants.hpp>
#include <stan/math/prim/fun/size.hpp>
#include <stan/math/prim/fun/size_zero.hpp>
#include <stan/math/prim/fun/value_of.hpp>
#include <stan/math/prim/functor/operands_and_partials.hpp>
namespace stan {
namespace math {
template <bool propto, typename T_y>
return_type_t<T_y> std_normal_lpdf(const T_y& y) {
static const char* function = "std_normal_lpdf";
using T_partials_return = partials_return_t<T_y>;
check_not_nan(function, "Random variable", y);
if (size_zero(y)) {
return 0.0;
}
if (!include_summand<propto, T_y>::value) {
return 0.0;
}
T_partials_return logp(0.0);
operands_and_partials<T_y> ops_partials(y);
scalar_seq_view<T_y> y_vec(y);
size_t N = stan::math::size(y);
for (size_t n = 0; n < N; n++) {
const T_partials_return y_val = value_of(y_vec[n]);
logp += y_val * y_val;
if (!is_constant_all<T_y>::value) {
ops_partials.edge1_.partials_[n] -= y_val;
}
}
logp *= -0.5;
if (include_summand<propto>::value) {
logp += NEG_LOG_SQRT_TWO_PI * N;
}
return ops_partials.build(logp);
}
template <typename T_y>
inline return_type_t<T_y> std_normal_lpdf(const T_y& y) {
return std_normal_lpdf<false>(y);
}
} // namespace math
} // namespace stan