Doc added

stan/math/prim/mat/prob/std_binormal_lcdf.hpp

@@ -1,23 +1,30 @@
-#ifndef STAN_MATH_PRIM_SCAL_PROB_STD_BINORMAL_LCDF_HPP
-#define STAN_MATH_PRIM_SCAL_PROB_STD_BINORMAL_LCDF_HPP
+#ifndef STAN_MATH_PRIM_MAT_PROB_STD_BINORMAL_LCDF_HPP
+#define STAN_MATH_PRIM_MAT_PROB_STD_BINORMAL_LCDF_HPP
 
 #include <stan/math/prim/scal/meta/is_constant_struct.hpp>
 #include <stan/math/prim/scal/meta/partials_return_type.hpp>
 #include <stan/math/prim/scal/meta/operands_and_partials.hpp>
 #include <stan/math/prim/scal/meta/scalar_seq_view.hpp>
+#include <stan/math/prim/mat/meta/operands_and_partials.hpp>
+#include <stan/math/prim/mat/meta/vector_seq_view.hpp>
 #include <stan/math/prim/scal/err/check_consistent_sizes.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/size_zero.hpp>
+#include <stan/math/prim/scal/err/check_size_match.hpp>
 #include <stan/math/prim/scal/fun/constants.hpp>
 #include <stan/math/prim/scal/fun/value_of.hpp>
+#include <stan/math/prim/scal/meta/size_of.hpp>
 #include <stan/math/prim/scal/meta/max_size.hpp>
 #include <stan/math/prim/scal/meta/contains_nonconstant_struct.hpp>
 #include <boost/random/normal_distribution.hpp>
 #include <stan/math/prim/scal/fun/Phi.hpp>
 #include <stan/math/prim/scal/prob/std_normal_lpdf.hpp>
 #include <stan/math/prim/scal/fun/std_binormal_integral.hpp>
 #include <stan/math/prim/scal/fun/inv.hpp>
+#include <stan/math/prim/scal/meta/max_size_mvt.hpp>
+#include <stan/math/prim/mat/meta/value_type.hpp>
+
 #include <cmath>
 #include <limits>
 
@@ -26,81 +33,98 @@ namespace math {
 
 /**
  * The log of the CDF of the standard bivariate normal (binormal) for the
- * specified scalar(s) given the specified correlation(s). y1, y2, or rho can
+ * specified vector(s) given the specified correlation(s). rho can
  * each be either a scalar or a vector. Any vector inputs must be the same
  * length.
  *
  * <p>The result log probability is defined to be the sum of the
- * log probabilities for each observation 1/observation 2/correlation 
- * parameter triple.
+ * log probabilities for each 2-vector observation/correlation 
+ * parameter pair.
  *
- * @param y1 (Sequence of) scalar(s).
- * @param y2 (Sequence of) scalar(s).
+ * @param y (Sequence of) 2-vector(s).
  * @param rho (Sequence of) correlation parameter(s)
  * for the standard bivariate normal distribution.
- * @return The log of the product of the probabilities P(Y1 <= y1, Y2 <= y2 | rho).
- * @throw std::domain_error if the rho is not between -1 and 1.
- * @tparam T_y1 Underlying type of random variable 1 in sequence.
- * @tparam T_y2 Underlying type of random variable 2 in sequence.
+ * @return The log of the product of the probabilities P(Y1 <= y[1], Y2 <= y[2] | rho).
+ * @throw std::domain_error if the rho is not between -1 and 1 or nan.
+ * @tparam T_y Underlying type of random variable in sequence.
  * @tparam T_rho Underlying type of correlation in sequence.
  */
 
-template <typename T_y1, typename T_y2, typename T_rho>
-typename return_type<T_y1, T_y2, T_rho>::type std_binormal_lcdf(
-    const T_y1& y1, const T_y2& y2, const T_rho& rho) {
+template <typename T_y, typename T_rho>
+typename return_type<T_y, T_rho>::type std_binormal_lcdf(
+    const T_y& y, const T_rho& rho) {
   static const char* function = "std_binormal_lcdf";
-  typedef typename stan::partials_return_type<T_y1, T_y2, T_rho>::type
+  typedef typename stan::partials_return_type<T_y, T_rho>::type
       T_partials_return;
+  typedef typename stan::math::value_type<T_y>::type T_y_child_type;
 
   using std::exp;
   using std::sqrt;
   using std::log;
   using std::asin;
+  using std::max;
+
+  check_bounded(function, "Correlation parameter", rho, -1.0, 1.0);
+  if (stan::is_vector_like<T_y_child_type>::value || stan::is_vector_like<T_rho>::value) {
+      check_consistent_sizes(function, "Random variable", y,
+                             "Correlation parameter", rho);
+  }
+
+
+  vector_seq_view<T_y> y_vec(y);
+  const scalar_seq_view<T_rho> rho_vec(rho);
+  size_t N_y = length_mvt(y);
+  size_t N_rho = size_of(rho);
+  size_t N = max(N_y, N_rho);
+
+  operands_and_partials<T_y, T_rho> ops_partials(y, rho);
 
   T_partials_return cdf_log(0.0);
-  if (size_zero(y1, y2, rho))
+  if (unlikely(N == 0))
     return cdf_log;
 
-  check_not_nan(function, "Random variable 1", y1);
-  check_not_nan(function, "Random variable 2", y2);
-  check_bounded(function, "Correlation parameter", rho, -1.0, 1.0);
-  check_consistent_sizes(function, "Random variable 1", y1, "Random variable 2",
-                         y2, "Correlation parameter", rho);
-
-  operands_and_partials<T_y1, T_y2, T_rho> ops_partials(y1, y2, rho);
+  // check size consistency of all random variables y
+  for (size_t i = 0; i < N_y; ++i) {
+    check_size_match(function,
+                     "Size of one of the vectors "
+                     "of the location variable",
+                     y_vec[i].size(),
+                     "Size of another vector of "
+                     "the location variable",
+                     2);
+  }
 
-  scalar_seq_view<T_y1> y1_vec(y1);
-  scalar_seq_view<T_y2> y2_vec(y2);
-  scalar_seq_view<T_rho> rho_vec(rho);
-  size_t N = max_size(y1, y2, rho);
+  for (size_t n = 0; n < N_y; ++n) {
+    check_not_nan(function, "Random variable", y_vec[n]);
+  }
 
   for (size_t n = 0; n < N; ++n) {
-    const T_partials_return y1_dbl = value_of(y1_vec[n]);
-    const T_partials_return y2_dbl = value_of(y2_vec[n]);
+    const T_partials_return y1_dbl = value_of(y_vec[n][0]);
+    const T_partials_return y2_dbl = value_of(y_vec[n][1]);
     const T_partials_return rho_dbl = value_of(rho_vec[n]);
 
     T_partials_return cdf_ = std_binormal_integral(y1_dbl, y2_dbl, rho_dbl);
       cdf_log += log(cdf_);
 
-    if (contains_nonconstant_struct<T_y1, T_y2, T_rho>::value) {
+    if (contains_nonconstant_struct<T_y, T_rho>::value) {
       const T_partials_return inv_cdf_ = cdf_ > 0 ? inv(cdf_) : 
         std::numeric_limits<double>::infinity();
       const T_partials_return one_minus_rho_sq = (1 + rho_dbl) * (1 - rho_dbl);
       const T_partials_return sqrt_one_minus_rho_sq = sqrt(one_minus_rho_sq);
       const T_partials_return rho_times_y2 = rho_dbl * y2_dbl;
       const T_partials_return y1_minus_rho_times_y2 = y1_dbl - rho_times_y2;
-      if (!is_constant_struct<T_y1>::value)
-        ops_partials.edge1_.partials_[n] 
+      if (!is_constant_struct<T_y>::value) {
+        ops_partials.edge1_.partials_vec_[n](0) 
           += cdf_ > 0 && cdf_ < 1 ? inv_cdf_ * exp(std_normal_lpdf(y1_dbl))
               * Phi((y2_dbl - rho_dbl * y1_dbl) / sqrt_one_minus_rho_sq)
               : cdf_ > 0 ? 1 : 0;
-      if (!is_constant_struct<T_y2>::value)
-        ops_partials.edge2_.partials_[n] 
+        ops_partials.edge1_.partials_vec_[n](1) 
           += cdf_ > 0 && cdf_ < 1 ? inv_cdf_ * exp(std_normal_lpdf(y2_dbl))
               * Phi(y1_minus_rho_times_y2 / sqrt_one_minus_rho_sq)
               : cdf_ > 0 ? 1 : 0;
+      }
       if (!is_constant_struct<T_rho>::value)
-        ops_partials.edge3_.partials_[n]
+        ops_partials.edge2_.partials_[n]
             += cdf_ > 0 && cdf_ < 1 ? inv_cdf_ * 0.5 / (stan::math::pi() * sqrt_one_minus_rho_sq)
             * exp(-0.5 / one_minus_rho_sq * y1_minus_rho_times_y2 *  y1_minus_rho_times_y2
                   -0.5 * y2_dbl * y2_dbl)

stan/math/prim/scal/fun/std_binormal_integral.hpp

@@ -21,6 +21,34 @@
 namespace stan {
 namespace math {
 
+/**
+ * The CDF of the standard bivariate normal (binormal) for the specified
+ * variable 1, z1, and variable 2, z2, given the specified correlation, rho.
+ *
+ * The function calculates the bivariate integral with Owen's algorithm for the
+ * domain fabs(z1) <= 1, fabs(z2) <= 1, and fabs(rho) < 0.7, and elsewhere uses
+ * the result std_binormal_integral(z1, z2, rho): 
+ *
+ * \f[
+ *   \Phi(z_1)\Phi(z_2) + 
+ *   \frac{1}{2\pi}\int_0^{\sin^{-1}(\rho)}\exp
+ *   (\tfrac{-z_1^2 - z_2^2 + 2z_1 z_2 \sin(\theta)}{2\cos(\theta)^2}) d\theta
+ * \f]
+ * 
+ * where the integral is done using boost's tanh-sinh quadrature. References
+ * for the integral can be found in:
+ *
+ * Genz, Alan. "Numerical computation of rectangular bivariate 
+ * and trivariate normal and t probabilities." 
+ * Statistics and Computing 14.3 (2004): 251-260.
+ *
+ * @param z1 Random variable 1
+ * @param z2 Random variable 2
+ * @param rho correlation parameter
+ * for the standard bivariate normal distribution.
+ * @return The probability P(Y1 <= y[1], Y2 <= y[2] | rho).
+ * @throw std::domain_error if the rho is not between -1 and 1 or nan.
+ */
 double std_binormal_integral(const double z1, 
                              const double z2, 
                              const double rho) {
@@ -76,9 +104,10 @@ double std_binormal_integral(const double z1,
                                                    &error, &L1);
     if (exp(log(L1) - log(accum)) > 1.5) {
       std::cout << "L1: " << L1 << " integral: " << accum << std::endl;
-      domain_error(function, "the numeric integration of bivariate normal density condition number", exp(log(L1) - log(accum)),
+      domain_error(function, "the numeric integration of bivariate normal density condition number", 
+                   exp(log(L1) - log(accum)),
                    "",
-                   " indicates the integral in poorly conditioned");
+                   " indicates the integral is poorly conditioned");
     }
 
 
@@ -91,6 +120,27 @@ double std_binormal_integral(const double z1,
   return z2 > -z1 > 0 ? Phi(z1) + Phi(z2) - 1 : 0;
 }
 
+/**
+ * The CDF of the standard bivariate normal (binormal) for the specified
+ * variable 1, z1, and variable 2, z2, given the specified correlation, rho.
+ * rho can each be either a scalar or a vector. If z1, z2, or rho are 
+ * autodiff types, the function propagates the gradients of the CDF. Reference
+ * for the gradient of the bivariate normal CDF can be found here:
+ *
+ * blogs.sas.com/content/iml/
+ * 2013/09/20/gradient-of-the-bivariate-normal-cumulative-distribution.html
+ *
+ * @param z1 Scalar random variable 1
+ * @param z2 Scalar random variable 2
+ * @param rho Scalar correlation parameter
+ * for the standard bivariate normal distribution.
+ * @return The probability P(Y1 <= y[1], Y2 <= y[2] | rho).
+ * @throw std::domain_error if the rho is not between -1 and 1 or nan.
+ * @tparam T_z1 type of z1
+ * @tparam T_z2 type of z2
+ * @tparam T_rho type of rho 
+ */
+
 template <typename T_z1, typename T_z2, typename T_rho>
 typename return_type<T_z1, T_z2, T_rho>::type
                                      std_binormal_integral(