Merge branch 'master' of https://code.google.com/p/stan

TO-DO.txt

@@ -15,25 +15,22 @@ The to-do list is organized into the following sections:
 V 1.0.2 Release Notes
 ======================================================================
 
-To Do for 1.0.2
-  ?? remove (unsigned) long long instances in dump and error testing
-  ?? manual section on reparameterization and change of vars
-
 Bug Fixes:
-  -- fixed sd() and variance() to return 0.0 for sequences of size 1
   -- check ranges for LHS of assignment to prevent seg faults
-  -- indexing fixed for arrays of matrix/vector/row_vector
-  -- vectorized several probability functions (see the manual)
-  -- fixed void return type in auto_covariance
-  -- removed template variable names from distribution error msgs
   -- added matrix size and shape tests to avoid seg faults
   -- changed matrix to throw domain_error rather than illegal_argument
+  -- removed template variable names from distribution error msgs
+  -- indexing fixed for arrays of matrix/vector/row_vector
+  -- fixed sd() and variance() to return 0.0 for sequences of size 1
+  -- fixed void return type in auto_covariance
   -- boundary conditions at limits of support for cdfs
   -- patch truncation to return -inf for variates out of range
   -- upgraded BUGS ring model to use constraints plus tan2()
 
 New Features:
   -- print statements
+  -- multiply_lower_tri_self_transpose function
+  -- vectorized several probability functions (see the manual)
 
 Manual Additions:
   -- programming guide: IRT models 
@@ -46,6 +43,8 @@ C++ API
 ======================================================================
 [items go here first, then to RStan/cmd]
 
+* remove (unsigned) long long instances in dump and error testing
+
 * bring isinf in line with revised Boost
 
 * add lower/upper bound infinite or negative infinite values

src/docs/stan-reference/acknowledgements.tex

@@ -41,9 +41,30 @@ \section*{Individuals}
 We thank Michael Betancourt for hanging out at our meetings and
 sharing his insights into differential geometry, John Salvatier for
 pointing us to automatic differentiation and \HMC in the first place,
-and Kristen van Leuven of {\sc iserp} for help preparing our grant
+and Kristen van Leuven of ISERP for help preparing our grant
 proposals.
 
+\subsection*{Code  and Doc Patches}
+
+Thanks to Jeffrey Arnold for submitting code patches and reporting
+a whole bunch of bugs for earlier versions. 
+
+Thanks to Jeffrey Arnold, Eric N.~Brown, and Wayne Folta for
+submitting documentation patches.
+
+\subsection*{Bug Reports}
+
+We're really thankful to everyone who's had the patience to try
+to get Stan working and reported bugs.  All the gory details are
+available from Stan's issue tracker at the following URL.
+%
+\begin{quote}
+\url{https://code.google.com/p/stan/issues/list}
+\end{quote}
+
+
+
+
 \vfill
 \begin{center}
 \hfill

src/docs/stan-reference/appendices.tex

@@ -1224,7 +1224,9 @@ \section{Line Length}
 
 Line lengths should not exceed 80 characters.%
 %
-\footnote{Even 80 characters may be too many for rendering in print.}
+\footnote{Even 80 characters may be too many for rendering in print;
+  for instance, in this manual, the number of code characters that fit
+  on a line is about 65.}
 %
 This is a typical recommendation for many programming language style
 guides because it makes it easier to lay out text edit windows side by
@@ -1256,7 +1258,7 @@ \section{Variable Naming}
 The exception to the lowercasing recommendation, which also follows
 the C/\Cpp conventions, is for size constants, for which the
 recommended form is a single uppercase letter.  The reason for this is
-that it allows the loop variables match.  So loops over the indices of
+that it allows the loop variables to match.  So loops over the indices of
 an $M \times N$ matrix $a$ would look as follows.
 %
 \begin{quote}

src/docs/stan-reference/commands.tex

@@ -900,7 +900,7 @@ \section{BNF Grammar for Dump Data}
 \begin{verbatim}
  definitions ::= definition+
 
- definition ::= name ("->" | '=') value optional_semicolon
+ definition ::= name ("<-" | '=') value optional_semicolon
 
  name ::= char* 
         | ''' char* ''' 
@@ -910,7 +910,7 @@ \section{BNF Grammar for Dump Data}
 
  value<T> ::= T 
             | seq<T>
-            | 'struct' '(' seq<T> ',' ".Dim" '=' seq<int> ')'
+            | 'structure' '(' seq<T> ',' ".Dim" '=' seq<int> ')'
 
  seq<int> ::= int ':' int
             | cseq<int>

src/docs/stan-reference/functions.tex

@@ -597,7 +597,7 @@ \subsection{Infix Matrix Operators}
 the scalar \farg{y} and vector \farg{x}}
 %
 \fitem{matrix}{operator*}{vector \farg{x}, row\_vector \farg{y}}{The product
-of the row vector \farg{x} and vector \farg{y}}
+of the vector \farg{x} and row vector \farg{y}}
 %
 %
 \fitem{row\_vector}{operator*}{row\_vector \farg{x}, real \farg{y}}{The product of

src/docs/stan-reference/programming.tex

@@ -443,9 +443,9 @@ \subsection{Estimating Censored Values}
 \end{Verbatim}
 \end{quote}
 %
-Because the censored data array \code{y} is declared to be a parameter, it
+Because the censored data array \code{y\_cens} is declared to be a parameter, it
 will be sampled along with the location and scale parameters \code{mu}
-and \code{sigma}.  Because the censored data array \code{y} is
+and \code{sigma}.  Because the censored data array \code{y\_cens} is
 declared to have values of type \code{real<lower=U>}, all imputed values
 for censored data will be greater than \code{U}.  The imputed censored
 data affects the location and scale parameters through the last
@@ -639,12 +639,13 @@ \section{Summing out the Responsibility Parameter}
 \end{quote}
 %
 The model involves \code{K} mixture components and \code{N} data
-points.  The mixing proportions are defined to be a unit $K$-simplex
-using \code{simplex[K]}, the components distributions locations
-\code{mu[k]} are unconstrained and their scales \code{sigma[k]}
-constrained to be positive.  The model declares a local variable
-\code{ps} of type \code{real[K]}, which is used to accumulate the
-contributions by each mixture component.
+points.  The mixing proportion parameter \code{theta} is declared to
+be a unit $K$-simplex, whereas the component location parameter
+\code{mu} and scale parameter \code{sigma} are both defined to be
+arrays of size \code{K}.  The values in the scale array \code{sigma}
+are constrained to be non-negative.  The model declares a local 
+array variable \code{ps} to be size \code{K} and uses it to
+accumulate the contributions from the mixture components.
 
 The locations and scales are drawn from simple priors for the sake of
 this example, but could be anything supported by \Stan.  The mixture
@@ -767,13 +768,13 @@ \section{Robust Noise Models}
 \begin{quote}
 \begin{Verbatim}
 data {
-    ...
-    real<lower=0> nu;
+  ...
+  real<lower=0> nu;
 }
 ...
 model {
-    for (n in 1:N)
-        y[n] ~ student_t(nu,0,sigma);
+  for (n in 1:N)
+    y[n] ~ student_t(nu, alpha + beta * x[n], sigma);
 }
 \end{Verbatim}
 \end{quote}
@@ -1158,13 +1159,13 @@ \subsection{Multilevel 2PL Model}
 \begin{quote}
 \begin{Verbatim}
 parameters {    
-  real delta;                    // mean student ability
-  real alpha[J];                 // ability for student j - mean ability
-  real beta[K];                  // difficulty for question k
-  real log_gamma[K];             // discriminativeness for question k
-  real<lower=0> sigma_alpha;     // sd of student abilities  
-  real<lower=0> sigma_beta;      // sd of question difficulties 
-  real<lower=0> sigma_gamma;     // sd of log question discriminativeness
+  real delta;                  // mean student ability
+  real alpha[J];               // ability for j - mean
+  real beta[K];                // difficulty for k
+  real log_gamma[K];           // discrim for k
+  real<lower=0> sigma_alpha;   // sd of abilities
+  real<lower=0> sigma_beta;    // sd of difficulties 
+  real<lower=0> sigma_gamma;   // sd of log discrim
 }
 \end{Verbatim}
 \end{quote}
@@ -1182,8 +1183,9 @@ \subsection{Multilevel 2PL Model}
   sigma_beta ~ cauchy(0,5);
   sigma_gamma ~ cauchy(0,5);
   for (n in 1:N)
-    y[n] ~ bernoulli_logit( exp(log_gamma[kk[n]])
-                            * (alpha[jj[n]] - beta[kk[n]] + delta) );
+    y[n] ~ bernoulli_logit( 
+               exp(log_gamma[kk[n]])
+               * (alpha[jj[n]] - beta[kk[n]] + delta) );
 }
 \end{Verbatim}
 \end{quote}

src/stan/agrad/matrix.hpp

@@ -1479,10 +1479,10 @@ namespace stan {
       int K = L.rows();
       matrix_v LLt(K,K);
       if (K == 0) return LLt;
-      if (K == 1) {
-        LLt(0,0) = L(0,0) * L(0,0);
-        return LLt;
-      }
+      // if (K == 1) {
+      //   LLt(0,0) = L(0,0) * L(0,0);
+      //   return LLt;
+      // }
       int Knz = (K * (K + 1)) / 2;  // nonzero: (K choose 2) below
                                     // diag + K on diag
       vari** vs = (vari**)memalloc_.alloc( Knz * sizeof(vari*) );

src/stan/command/stanc.cpp

@@ -155,7 +155,7 @@ int main(int argc, const char* argv[]) {
     std::cout << "Output file=" << out_file_name << std::endl;
     bool include_main = !cmd.has_flag("no_main");
     bool valid_model 
-      = stan::gm::compile(in,out,model_name,include_main);
+      = stan::gm::compile(in,out,model_name,include_main,in_file_name);
     out.close();
     if (!valid_model) {
       std::cout << "PARSING FAILED." << std::endl;

src/stan/gm/compiler.hpp

@@ -20,18 +20,21 @@ namespace stan {
      * @param stan_gm_in Graphical model specification
      * @param cpp_out C++ code output stream
      * @param model_name Name of model class
-     * @param include_main Indicates whether to generate a 
-     *    main function
+     * @param include_main Indicates whether to generate a main
+     *    function
+     * @param in_file_name Name of input file to use in error
+     * messages; defaults to <code>input</code>.
      * @return <code>false</code> if code could not be generated
      *    due to syntax error in the Graphical model; 
      *    <code>true</code> otherwise.
      */
     bool compile(std::istream& stan_gm_in,
                  std::ostream& cpp_out,
                  const std::string& model_name,
-                 bool include_main = true) {
+                 bool include_main = true,
+                 const std::string& in_file_name = "input") {
       program prog;
-      bool parsed_ok = parse(stan_gm_in,"input",prog);
+      bool parsed_ok = parse(stan_gm_in,in_file_name,prog);
       if (!parsed_ok) 
         return false; // syntax error in program
       generate_cpp(prog,model_name,cpp_out,include_main);

src/stan/prob/distributions/univariate/continuous/beta.hpp

@@ -40,7 +40,8 @@ namespace stan {
     beta_log(const T_y& y, const T_scale_succ& alpha, const T_scale_fail& beta, 
              const Policy&) {
       static const char* function = "stan::prob::beta_log(%1%)";
-
+
+      using stan::is_constant_struct;
       using stan::math::check_positive;
       using stan::math::check_finite;
       using stan::math::check_not_nan;
@@ -104,49 +105,51 @@ namespace stan {
       agrad::OperandsAndPartials<T_y, T_scale_succ, T_scale_fail> operands_and_partials(y, alpha, beta);
 
 
-      DoubleVectorView<true,T_y> log_y(length(y));
-      for (size_t n = 0; n < length(y); n++)
-	log_y[n] = log(value_of(y_vec[n]));
+      DoubleVectorView<include_summand<propto,T_y,T_scale_succ>::value,T_y> log_y(length(y));
       DoubleVectorView<include_summand<propto,T_y,T_scale_fail>::value,T_y> log1m_y(length(y));
-      if (include_summand<propto,T_y,T_scale_fail>::value) {
-	for (size_t n = 0; n < length(y); n++)
+
+      for (size_t n = 0; n < length(y); n++) {
+	if (include_summand<propto,T_y,T_scale_succ>::value)
+	  log_y[n] = log(value_of(y_vec[n]));
+	if (include_summand<propto,T_y,T_scale_fail>::value)
 	  log1m_y[n] = log1m(value_of(y_vec[n]));
       }
 
-      DoubleVectorView<include_summand<propto,T_scale_succ>::value,T_scale_succ> lgamma_alpha(length(alpha));
-      DoubleVectorView<include_summand<propto,T_scale_succ>::value,T_scale_succ> digamma_alpha(length(alpha));
-      if (include_summand<propto,T_scale_succ>::value) {
-	for (size_t n = 0; n < length(alpha); n++) {
+    DoubleVectorView<include_summand<propto,T_scale_succ>::value,T_scale_succ> lgamma_alpha(length(alpha));
+      DoubleVectorView<!is_constant_struct<T_scale_succ>::value,T_scale_succ> digamma_alpha(length(alpha));
+      for (size_t n = 0; n < length(alpha); n++) {
+	if (include_summand<propto,T_scale_succ>::value) 
 	  lgamma_alpha[n] = lgamma(value_of(alpha_vec[n]));
+	if (!is_constant_struct<T_scale_succ>::value)
 	  digamma_alpha[n] = digamma(value_of(alpha_vec[n]));
-	}
       }
 
+
       DoubleVectorView<include_summand<propto,T_scale_fail>::value,T_scale_fail> lgamma_beta(length(beta));
-      DoubleVectorView<include_summand<propto,T_scale_fail>::value,T_scale_fail> digamma_beta(length(beta));
-      if (include_summand<propto,T_scale_fail>::value) {
-	for (size_t n = 0; n < length(beta); n++) {
+      DoubleVectorView<!is_constant_struct<T_scale_fail>::value,T_scale_fail> digamma_beta(length(beta));
+      for (size_t n = 0; n < length(beta); n++) {
+	if (include_summand<propto,T_scale_fail>::value) 
 	  lgamma_beta[n] = lgamma(value_of(beta_vec[n]));
+	if (!is_constant_struct<T_scale_fail>::value)
 	  digamma_beta[n] = digamma(value_of(beta_vec[n]));
-	}
       }
 
       DoubleVectorView<include_summand<propto,T_scale_succ,T_scale_fail>::value,
 	double,
 	is_vector<T_scale_succ>::value||is_vector<T_scale_fail>::value>
 	lgamma_alpha_beta(max_size(alpha,beta));
-      DoubleVectorView<include_summand<propto,T_scale_succ,T_scale_fail>::value,
+      DoubleVectorView<!is_constant_struct<T_scale_succ>::value||!is_constant_struct<T_scale_fail>::value,
 	double,
 	is_vector<T_scale_succ>::value||is_vector<T_scale_fail>::value>
 	digamma_alpha_beta(max_size(alpha,beta));
-      if (include_summand<propto,T_scale_succ,T_scale_fail>::value) {
-	for (size_t n = 0; n < max_size(alpha,beta); n++) {
-	  lgamma_alpha_beta[n] = lgamma(value_of(alpha_vec[n]) + value_of(beta_vec[n]));
-	  digamma_alpha_beta[n] = digamma(value_of(alpha_vec[n]) + value_of(beta_vec[n]));
-	}
+      for (size_t n = 0; n < max_size(alpha,beta); n++) {
+	const double alpha_beta = value_of(alpha_vec[n]) + value_of(beta_vec[n]);
+	if (include_summand<propto,T_scale_succ,T_scale_fail>::value)
+	  lgamma_alpha_beta[n] = lgamma(alpha_beta);
+	if (!is_constant_struct<T_scale_succ>::value||!is_constant_struct<T_scale_fail>::value)
+	  digamma_alpha_beta[n] = digamma(alpha_beta);
       }
 
-
       for (size_t n = 0; n < N; n++) {
 	// pull out values of arguments
 	const double y_dbl = value_of(y_vec[n]);
@@ -166,11 +169,11 @@ namespace stan {
 	  logp += (beta_dbl-1.0) * log1m_y[n];
 
 	// gradients
-	if (include_summand<propto,T_y>::value)
+	if (!is_constant_struct<T_y>::value)
 	  operands_and_partials.d_x1[n] += (alpha_dbl-1)/y_dbl + (beta_dbl-1)/(y_dbl-1);
-	if (include_summand<propto,T_scale_succ>::value)
+	if (!is_constant_struct<T_scale_succ>::value)
 	  operands_and_partials.d_x2[n] += log_y[n] + digamma_alpha_beta[n] - digamma_alpha[n];
-	if (include_summand<propto,T_scale_fail>::value)
+	if (!is_constant_struct<T_scale_fail>::value)
 	  operands_and_partials.d_x3[n] += log1m_y[n] + digamma_alpha_beta[n] - digamma_beta[n];
       }
       return operands_and_partials.to_var(logp);

src/stan/prob/distributions/univariate/continuous/cauchy.hpp

@@ -36,7 +36,8 @@ namespace stan {
     cauchy_log(const T_y& y, const T_loc& mu, const T_scale& sigma, 
                const Policy&) {
       static const char* function = "stan::prob::cauchy_log(%1%)";
-
+
+      using stan::is_constant_struct;
       using stan::math::check_positive;
       using stan::math::check_finite;
       using stan::math::check_not_nan;
@@ -121,11 +122,11 @@ namespace stan {
 	  logp -= log1p(y_minus_mu_over_sigma_squared);
 
         // gradients
-	if (!is_constant<typename is_vector<T_y>::type>::value)
+	if (!is_constant_struct<T_y>::value)
 	  operands_and_partials.d_x1[n] -= 2 * y_minus_mu / (sigma_squared[n] + y_minus_mu_squared);
-        if (!is_constant<typename is_vector<T_loc>::type>::value)
+        if (!is_constant_struct<T_loc>::value)
           operands_and_partials.d_x2[n] += 2 * y_minus_mu / (sigma_squared[n] + y_minus_mu_squared);
-        if (!is_constant<typename is_vector<T_scale>::type>::value)
+        if (!is_constant_struct<T_scale>::value)
           operands_and_partials.d_x3[n] += (y_minus_mu_squared - sigma_squared[n]) * inv_sigma[n] / (sigma_squared[n] + y_minus_mu_squared);
       }
       return operands_and_partials.to_var(logp);