Thanks @Foztarz. I ended up implementing a C++ version of unit vector von Mises lpdf myself. I am using it with rstan. Here is the code:

```
#include <stan/model/model_base.hpp>
#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/scal/err/check_consistent_sizes.hpp>
#include <stan/math/prim/scal/err/check_finite.hpp>
#include <stan/math/prim/scal/err/check_positive_finite.hpp>
#include <stan/math/prim/scal/fun/log_modified_bessel_first_kind.hpp>
#include <stan/math/prim/scal/fun/modified_bessel_first_kind.hpp>
#include <stan/math/prim/scal/fun/constants.hpp>
#include <stan/math/prim/scal/fun/value_of.hpp>
#include <stan/math/prim/scal/fun/size_zero.hpp>
#include <cmath>
#include <boost/math/tools/promotion.hpp>
using namespace stan;
template <bool propto, typename T_y, typename T_mu, typename T_kappa>
typename boost::math::tools::promote_args<T_y, T_mu, T_kappa>::type
unit_vector_von_mises_lpdf(
const Eigen::Matrix<T_y, Eigen::Dynamic, Eigen::Dynamic>& y,
const Eigen::Matrix<T_mu, 1, Eigen::Dynamic>& mu,
const Eigen::Matrix<T_kappa, Eigen::Dynamic, 1>& kappa,
std::ostream* pstream__) {
if (y.rows() == 0) {
return 0;
}
typedef typename stan::partials_return_type<
Eigen::Matrix<T_y, Eigen::Dynamic, Eigen::Dynamic>,
Eigen::Matrix<T_mu, 1, Eigen::Dynamic>,
Eigen::Matrix<T_kappa, Eigen::Dynamic, 1>>::type T_partials_return;
T_partials_return logp(0);
const bool y_const = is_constant_all<
Eigen::Matrix<T_y, Eigen::Dynamic, Eigen::Dynamic>>::value;
const bool mu_const = is_constant_all<
Eigen::Matrix<T_mu, 1, Eigen::Dynamic>>::value;
const bool kappa_const = is_constant_all<
Eigen::Matrix<T_kappa, Eigen::Dynamic, 1>>::value;
operands_and_partials<
Eigen::Matrix<T_y, Eigen::Dynamic, Eigen::Dynamic>,
Eigen::Matrix<T_mu, 1, Eigen::Dynamic>,
Eigen::Matrix<T_kappa, Eigen::Dynamic, 1>> ops_partials(y, mu, kappa);
for (size_t i = 0; i < y.rows(); i++) {
Eigen::Matrix<T_y, Eigen::Dynamic, 1> y_row(y.row(i));
const T_partials_return kappa_dbl = value_of(kappa(i));
const T_partials_return mu_dot_y = value_of((mu*y_row)(0));
if (include_summand<propto>::value) {
logp -= LOG_TWO_PI;
}
if (include_summand<propto, T_kappa>::value) {
logp -= log_modified_bessel_first_kind(0, kappa_dbl);
}
if (include_summand<propto, T_y, T_mu, T_kappa>::value) {
logp += kappa_dbl*mu_dot_y;
}
if (!y_const) {
for (size_t j = 0; j < length(mu); j++) {
ops_partials.edge1_.partials_vec_[i][j] += kappa_dbl*value_of(mu(j));
}
}
if (!mu_const) {
for (size_t j = 0; j < length(y_row); j++) {
ops_partials.edge2_.partials_vec_[i][j] += kappa_dbl*value_of(y_row(j));
}
}
if (!kappa_const) {
const T_partials_return bessel1 =
modified_bessel_first_kind(-1, kappa_dbl);
const T_partials_return bessel0 =
modified_bessel_first_kind(0, kappa_dbl);
ops_partials.edge3_.partials_[i] += mu_dot_y - bessel1/bessel0;
}
}
return ops_partials.build(logp);
}
template <bool propto, typename T_y, typename T_mu, typename T_kappa>
typename boost::math::tools::promote_args<T_y, T_mu, T_kappa>::type
unit_vector_von_mises_lpdf(
const Eigen::Matrix<T_y, 1, Eigen::Dynamic>& y,
const Eigen::Matrix<T_mu, 1, Eigen::Dynamic>& mu,
const T_kappa& kappa,
std::ostream* pstream__) {
using T_partials_return = partials_return_t<
Eigen::Matrix<T_y, 1, Eigen::Dynamic>,
Eigen::Matrix<T_mu, 1, Eigen::Dynamic>,
T_kappa>;
if (size_zero(y, mu, kappa)) {
return 0.0;
}
using std::log;
T_partials_return logp = 0.0;
if (!include_summand<propto, T_y, T_mu, T_kappa>::value) {
return logp;
}
const bool y_const = is_constant_all<
Eigen::Matrix<T_y, 1, Eigen::Dynamic>>::value;
const bool mu_const = is_constant_all<
Eigen::Matrix<T_mu, 1, Eigen::Dynamic>>::value;
const bool kappa_const = is_constant_all<T_kappa>::value;
operands_and_partials<
Eigen::Matrix<T_y, 1, Eigen::Dynamic>,
Eigen::Matrix<T_mu, 1, Eigen::Dynamic>,
T_kappa> ops_partials(y, mu, kappa);
Eigen::Matrix<T_y, Eigen::Dynamic, 1> y_row(y.transpose());
const T_partials_return kappa_dbl = value_of(kappa);
const T_partials_return mu_dot_y = value_of((mu*y_row)(0));
if (include_summand<propto>::value) {
logp -= LOG_TWO_PI;
}
if (include_summand<propto, T_kappa>::value) {
logp -= log_modified_bessel_first_kind(0, kappa_dbl);
}
if (include_summand<propto, T_y, T_mu, T_kappa>::value) {
logp += kappa_dbl*mu_dot_y;
}
if (!y_const) {
for (size_t j = 0; j < length(y_row); j++) {
ops_partials.edge1_.partials_[j] += kappa_dbl*value_of(mu(j));
}
}
if (!mu_const) {
for (size_t j = 0; j < length(mu); j++) {
ops_partials.edge2_.partials_[j] += kappa_dbl*value_of(y_row(j));
}
}
if (!kappa_const) {
const T_partials_return bessel1 = modified_bessel_first_kind(-1, kappa_dbl);
const T_partials_return bessel0 = modified_bessel_first_kind(0, kappa_dbl);
ops_partials.edge3_.partials_[0] += mu_dot_y - bessel1/bessel0;
}
return ops_partials.build(logp);
}
template <typename T_y, typename T_mu, typename T_kappa>
inline typename boost::math::tools::promote_args<T_y, T_mu, T_kappa>::type
unit_vector_von_mises_lpdf(const Eigen::Matrix<T_y, 1, Eigen::Dynamic>& y,
const Eigen::Matrix<T_mu, 1, Eigen::Dynamic>& mu,
const T_kappa& kappa, std::ostream* pstream__) {
return unit_vector_von_mises_lpdf<false>(y, mu, kappa, pstream__);
}
```