Commit 3cadc9dd authored by Ben Goodrich's avatar Ben Goodrich
Browse files

introduce inv_Phi()

parent 01c2a912
No related merge requests found
Showing with 331 additions and 1 deletion
+331 -1
# /
/demo
/models
/test
/*.dSYM
/bin
......
#ifndef STAN_MATH_FWD_SCAL_FUN_INV_PHI_HPP
#define STAN_MATH_FWD_SCAL_FUN_INV_PHI_HPP
#include <stan/math/fwd/core.hpp>
#include <stan/math/prim/scal/fun/inv_Phi.hpp>
#include <stan/math/prim/scal/fun/constants.hpp>
namespace stan {
namespace math {
template <typename T>
inline fvar<T> inv_Phi(const fvar<T>& x) {
using stan::math::inv_Phi;
using std::exp;
using std::log;
T xv = inv_Phi(x.val_);
return fvar<T>(xv,
x.d_ / exp(-0.5 * xv * xv) * SQRT_2_TIMES_SQRT_PI);
}
}
}
#endif
#ifndef STAN_MATH_PRIM_MAT_FUN_INV_PHI_HPP
#define STAN_MATH_PRIM_MAT_FUN_INV_PHI_HPP
#include <stan/math/prim/mat/fun/Eigen.hpp>
#include <stan/math/prim/mat/fun/inv_Phi.hpp>
namespace stan {
namespace math {
/**
* Return the element-wise standard normal inverse CDF
*
* @param p The vector / matrix of cumulative probabilities
* @return ret(i, j) = inv_Phi(m(i, j))
*/
template<typename T, int Rows, int Cols>
inline Eigen::Matrix<T, Rows, Cols>
inv_Phi(const Eigen::Matrix<T, Rows, Cols>& m) {
using stan::math::inv_Phi;
Eigen::Matrix<T, Rows, Cols> ret(m.rows(), m.cols());
for (size_t j = 0; j < m.cols(); j++)
for (size_t i = 0; i < m.rows(); i++)
ret.coeffRef(i,j) = inv_Phi(m.coeffRef(i,j));
return ret;
}
}
}
#endif
#ifndef STAN_MATH_PRIM_SCAL_FUN_INV_PHI_HPP
#define STAN_MATH_PRIM_SCAL_FUN_INV_PHI_HPP
#include <boost/math/tools/promotion.hpp>
#include <boost/math/special_functions/erf.hpp>
#include <stan/math/prim/scal/fun/constants.hpp>
#include <stan/math/prim/scal/err/check_not_nan.hpp>
#include <stan/math/prim/scal/err/check_bounded.hpp>
#include <stan/math/prim/scal/fun/Phi.hpp>
#include <stan/math/prim/scal/fun/log1m.hpp>
namespace stan {
namespace math {
/**
* The inverse of the unit normal cumulative distribution function.
*
* The return value for a specified input probability, $p$, is
* the unit normal variate, $x$, such that
*
* \f$\Phi(x) = \int_{-\infty}^x \mbox{\sf Norm}(x|0, 1) \ dx = p\f$
*
* Algorithm first derived in 2003 by Peter Jon Aklam at
* http://home.online.no/~pjacklam/notes/invnorm/
*
* @param p Argument between 0 and 1.
* @return Real number
*/
template <typename T>
inline typename boost::math::tools::promote_args<T>::type
inv_Phi(const T p) {
// overridden in fvar and var, so can hard-code
// here for scalars only
using stan::math::check_not_nan;
check_not_nan("inv_Phi", "p", p);
stan::math::check_bounded<T, double, double>
("inv_Phi", "Probability variable", p, 0, 1);
// if T is not double, the gradient will overflow before inv_Phi does
if (p < 8e-311)
return NEGATIVE_INFTY;
if (p == 1)
return INFTY;
const double a[6] = {
-3.969683028665376e+01, 2.209460984245205e+02,
-2.759285104469687e+02, 1.383577518672690e+02,
-3.066479806614716e+01, 2.506628277459239e+00
};
const double b[5] = {
-5.447609879822406e+01, 1.615858368580409e+02,
-1.556989798598866e+02, 6.680131188771972e+01,
-1.328068155288572e+01
};
const double c[6] = {
-7.784894002430293e-03, -3.223964580411365e-01,
-2.400758277161838e+00, -2.549732539343734e+00,
4.374664141464968e+00, 2.938163982698783e+00
};
const double d[4] = {
7.784695709041462e-03, 3.224671290700398e-01,
2.445134137142996e+00, 3.754408661907416e+00
};
const double p_low = 0.02425;
const double p_high = 0.97575;
T q, r, x;
if ( (p_low <= p) && (p <= p_high) ) { // central region
q = p - 0.5;
r = q * q;
x = (((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r+a[5])*q /
(((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r+1.0);
} else if (p < p_low) { // lower region
q = std::sqrt(-2.0*std::log(p));
x = (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) /
((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1.0);
} else { // upper region
q = std::sqrt(-2.0*stan::math::log1m(p));
x = -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) /
((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1.0);
}
if (x < 37.6) { // gradient blows up past here anyway
T e, u;
e = stan::math::Phi(x) - p;
u = e * SQRT_2_TIMES_SQRT_PI * std::exp(0.5 * x * x);
x = x - u / (1.0 + 0.5 * x * u);
}
return x;
}
}
}
#endif
#ifndef STAN_MATH_REV_SCAL_FUN_INV_PHI_HPP
#define STAN_MATH_REV_SCAL_FUN_INV_PHI_HPP
#include <stan/math/rev/core.hpp>
#include <stan/math/prim/scal/fun/inv_Phi.hpp>
#include <stan/math/prim/scal/fun/constants.hpp>
namespace stan {
namespace math {
namespace {
class inv_Phi_vari : public op_v_vari {
private:
double m_z;
public:
explicit inv_Phi_vari(vari* avi) :
op_v_vari(m_z = stan::math::inv_Phi(avi->val_), avi) {
}
void chain() {
static const double NEG_HALF = -0.5;
avi_->adj_ += adj_
* SQRT_2_TIMES_SQRT_PI
/ std::exp(NEG_HALF * m_z * m_z);
}
};
}
/**
* The inverse of unit normal cumulative density function.
*
* See stan::math::inv_Phi() for the double-based version.
*
* The derivative is the reciprocal of unit normal density function,
*
* @param p Probability
* @return The unit normal inverse cdf evaluated at p
*/
inline var inv_Phi(const stan::math::var& a) {
return var(new inv_Phi_vari(a.vi_));
}
}
}
#endif
#include <gtest/gtest.h>
#include <stan/math/fwd/scal/fun/inv_Phi.hpp>
#include <stan/math/prim/scal/fun/inv_Phi.hpp>
#include <stan/math/prim/scal/prob/normal_log.hpp>
#include <test/unit/math/fwd/scal/fun/nan_util.hpp>
#include <stan/math/fwd/core.hpp>
#include <stan/math/fwd/scal/fun/exp.hpp>
TEST(AgradFwdinv_Phi,Fvar) {
using stan::math::fvar;
using stan::math::inv_Phi;
fvar<double> x = 0.1;
x.d_ = 1.0;
fvar<double> inv_Phi_x = inv_Phi(x);
EXPECT_FLOAT_EQ(inv_Phi(0.1), inv_Phi_x.val_);
EXPECT_FLOAT_EQ(1.0 / exp(stan::math::normal_log<false>(inv_Phi_x.val_,0.0,1.0)),
inv_Phi_x.d_);
}
TEST(AgradFwdinv_Phi, FvarFvarDouble) {
using stan::math::fvar;
using stan::math::inv_Phi;
fvar<fvar<double> > x;
x.val_.val_ = 0.1;
x.val_.d_ = 1.0;
fvar<fvar<double> > a = inv_Phi(x);
EXPECT_FLOAT_EQ(inv_Phi(0.1), a.val_.val_);
EXPECT_FLOAT_EQ(1.0 / exp(stan::math::normal_log<false>(a.val_.val_,0.0,1.0)), a.val_.d_);
EXPECT_FLOAT_EQ(0, a.d_.val_);
EXPECT_FLOAT_EQ(0, a.d_.d_);
fvar<fvar<double> > y;
y.val_.val_ = 0.1;
y.d_.val_ = 1.0;
a = inv_Phi(y);
EXPECT_FLOAT_EQ(inv_Phi(0.1), a.val_.val_);
EXPECT_FLOAT_EQ(0, a.val_.d_);
EXPECT_FLOAT_EQ(1.0 / exp(stan::math::normal_log<false>(a.val_.val_,0.0,1.0)), a.d_.val_);
EXPECT_FLOAT_EQ(0, a.d_.d_);
}
struct inv_Phi_fun {
template <typename T0>
inline T0
operator()(const T0& arg1) const {
return inv_Phi(arg1);
}
};
TEST(AgradFwdinv_Phi,inv_Phi_NaN) {
inv_Phi_fun inv_Phi_;
test_nan_fwd(inv_Phi_,true);
}
#include <boost/math/special_functions/fpclassify.hpp>
#include <stan/math/prim/scal/fun/inv_Phi.hpp>
#include <stan/math/prim/scal/fun/Phi.hpp>
#include <gtest/gtest.h>
TEST(MathFunctions, inv_Phi) {
using stan::math::inv_Phi;
using stan::math::Phi;
EXPECT_FLOAT_EQ(0.0, inv_Phi(0.5));
double p = 0.123456789;
EXPECT_FLOAT_EQ(p, Phi(inv_Phi(p)));
p = 8e-311;
EXPECT_FLOAT_EQ(p, Phi(inv_Phi(p)));
p = 0.99;
EXPECT_FLOAT_EQ(p, Phi(inv_Phi(p)));
}
TEST(MathFunctions, Phi_inf) {
using stan::math::inv_Phi;
double p = 7e-311;
const double inf = std::numeric_limits<double>::infinity();
EXPECT_EQ(inv_Phi(p),-inf);
p = 0.99999999999999999999999999999999999999;
EXPECT_EQ(inv_Phi(p),inf);
}
TEST(MathFunctions, Phi_nan) {
using stan::math::inv_Phi;
double nan = std::numeric_limits<double>::quiet_NaN();
EXPECT_THROW(inv_Phi(nan), std::domain_error);
EXPECT_THROW(inv_Phi(-2.0), std::domain_error);
EXPECT_THROW(inv_Phi(2.0), std::domain_error);
}
#include <stan/math/rev/scal/fun/inv_Phi.hpp>
#include <test/unit/math/rev/mat/fun/util.hpp>
#include <gtest/gtest.h>
#include <test/unit/math/rev/scal/fun/nan_util.hpp>
TEST(AgradRev, inv_Phi) {
using stan::math::var;
std::vector<double> p_values;
p_values.push_back(0.1);
p_values.push_back(0.5);
p_values.push_back(0.75);
// d/dp = sqrt(2 * pi) / exp(normal_log(value_of(x), 0.0, 1.0))
std::vector<double> dp_values;
dp_values.push_back(5.69805985611701);
dp_values.push_back(2.506628274631);
dp_values.push_back(3.14686508056109);
for (size_t i = 0; i < p_values.size(); i++) {
var p, inv_Phi_p;
AVEC x;
VEC dp;
p = p_values[i];
inv_Phi_p = stan::math::inv_Phi(p);
x = createAVEC(p);
inv_Phi_p.grad(x,dp);
EXPECT_FLOAT_EQ(stan::math::inv_Phi(p.val()), inv_Phi_p.val());
EXPECT_FLOAT_EQ(dp_values[i], dp[0])
<< "p = " << p;
}
}
struct inv_Phi_fun {
template <typename T0>
inline T0
operator()(const T0& arg1) const {
return inv_Phi(arg1);
}
};
TEST(AgradRev,inv_Phi_NaN) {
inv_Phi_fun inv_Phi_;
test_nan(inv_Phi_,true,false);
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment