stan-dev / math Goto Github PK

From @betanalpha on August 16, 2013 14:28

Expose Eigen's FFT into the modeling language:

vector fft(vector x)

Copied from original issue: stan-dev/stan#179

List of functions:

• stan/math/prim/arr/err/check_matching_sizes.hpp
• stan/math/prim/arr/err/check_nonzero_size.hpp
• stan/math/prim/arr/err/check_ordered.hpp
• stan/math/prim/mat/err/check_cholesky_factor.hpp
• stan/math/prim/mat/err/check_cholesky_factor_corr.hpp
• stan/math/prim/mat/err/check_column_index.hpp
• stan/math/prim/mat/err/check_corr_matrix.hpp
• stan/math/prim/mat/err/check_cov_matrix.hpp
• stan/math/prim/mat/err/check_ldlt_factor.hpp
• stan/math/prim/mat/err/check_lower_triangular.hpp
• stan/math/prim/mat/err/check_matching_dims.hpp
• stan/math/prim/mat/err/check_multiplicable.hpp
• stan/math/prim/mat/err/check_ordered.hpp
• stan/math/prim/mat/err/check_pos_definite.hpp
• stan/math/prim/mat/err/check_pos_semidefinite.hpp
• stan/math/prim/mat/err/check_positive_ordered.hpp
• stan/math/prim/mat/err/check_range.hpp
• stan/math/prim/mat/err/check_row_index.hpp
• stan/math/prim/mat/err/check_simplex.hpp
• stan/math/prim/mat/err/check_spsd_matrix.hpp
• stan/math/prim/mat/err/check_square.hpp
• stan/math/prim/mat/err/check_std_vector_index.hpp
• stan/math/prim/mat/err/check_symmetric.hpp
• stan/math/prim/mat/err/check_unit_vector.hpp
• stan/math/prim/mat/err/check_vector.hpp
• stan/math/prim/scal/err/check_bounded.hpp
• stan/math/prim/scal/err/check_finite.hpp
• stan/math/prim/scal/err/check_greater.hpp
• stan/math/prim/scal/err/check_greater_or_equal.hpp
• stan/math/prim/scal/err/check_less.hpp
• stan/math/prim/scal/err/check_less_or_equal.hpp
• stan/math/prim/scal/err/check_nonnegative.hpp
• stan/math/prim/scal/err/check_not_nan.hpp
• stan/math/prim/scal/err/check_positive.hpp
• stan/math/prim/scal/err/check_positive_finite.hpp
• stan/math/prim/scal/err/check_positive_size.hpp
• stan/math/prim/scal/err/check_size_match.hpp

Not applicable

• stan/math/prim/scal/err/check_consistent_size.hpp
• stan/math/prim/scal/err/check_consistent_sizes.hpp

From @syclik on October 26, 2014 6:44

We want simple functions that can be exposed through the language.

Copied from original issue: stan-dev/stan#1104

From @syclik on November 22, 2013 16:21

Issue #333 shows the gradregIncGamma function getting into an infinite loop. The additional functions need to be tested:

• F32
• gradF32
• grad2F1
• gradRegIncBeta

Copied from original issue: stan-dev/stan#410

From @bob-carpenter on September 19, 2014 21:57

It would be more efficient to have densities that can be used like

cholesky_factor_cov[K] L;
L ~ wishart_cholesky(..,L0);
Sigma <- L * L';

Where Sigma would then have a Wishart density and L0 would be the cholesky factor of the parameter.

Copied from original issue: stan-dev/stan#1037

From @bob-carpenter on August 27, 2014 17:1

From Ben Goodrich on stan-dev:

It is kind of like the ODE thing but going the other way. In the case of a copula, we know the explicit form of the multivariate CDF of the marginal uniforms

F(u_1, u_2, ..., u_d | theta)

but for Stan's purposes we need the (logarithm of the) copula density function, which is

dF(u_1, u_2, ..., u_d | theta)
------------------------------ = f(u_1, u_2, ..., u_d | theta)
du_1,du_2,...,du_d

but for large d, there is sometimes no explicit form of f(u_1, u_2, ..., u_d | theta). So, loosely speaking, it would be great if we could evaluate

var f_log(var_vector u, var theta) {
  const int d = u.rows();
  var p = F(u, theta);
  var dp_du_1 = get_partial(p, u[1]);
  var dp_du_1_du_2 = get_partial(dp_du_1, u[2]);
  /* more of these */
  var log_density = log(get_partial(dp_du_*, u[d]);
  return log_density;
}

I replied:

We can do this with second-order autodiff. But it's not going to be super efficient in high dimensions, because each dimension is going to involve a forward pass on the function F, unless we can figure out how to do nested evaluation the way we are building up the ode solver.

The get_partial() operation you have will be a forward-mode auto-diff d_ value. That can then be logged, or whatever, and autodiffed through for HMC.

If the F doesn't involve any pdf or cdf calls in Stan, could do it now. Otherwise, we're going to have to wait until the 2nd order autodiff is completed. It's held up on Windows testing at the moment.

Are there a fixed number of these, or do people write their own F?

To which Ben replied:

There are potentially an infinite number of copula CDFs but in practice, people can only use the ones that have an explicit copula density function for estimation. But also in practice, they estimate the parameters of each margin separately and then pretend those estimates are fixed when they estimate the copula part. Full Bayesian estimation would be much better if we could make it feasible in high dimensions.

Copied from original issue: stan-dev/stan#938

From @bob-carpenter on January 11, 2015 22:39

Suggested by Joshua N Pritikin via stan-dev:

The traits metaprograms stan::math::value_type and stan::scalar_type are near duplicates and should be merged.

• The code layout for stan::math::value_type is the right one. It goes in math and separates matrices from standard library uses. It probably has the right test setup, too, since I wrote it recently.
• The extension to pointers in stan::scalar_type is going to have to go in the combined definition.

I'm not sure what all the helper-class complexity is for in stan::scalar_type but it'd be nice to understand before changing the behavior.

Copied from original issue: stan-dev/stan#1211

From @bob-carpenter on January 3, 2015 2:22

Need utilities for Kroneckers for Gaussian processes. I don't know much about it but want to get it down as an issue and point to Ben's R example:

https://groups.google.com/forum/#!msg/stan-users/iddShgJVwTQ/QgcDEP3h9JQJ

Talk to Ben to figure out how this should look in Stan.

Copied from original issue: stan-dev/stan#1196

From @chjackson on June 8, 2014 13:47

Feature request for a matrix exponential function, as raised recently on the stan-users list.

Copied from original issue: stan-dev/stan#680

From @jpritikin on January 5, 2015 21:10

This compiles but explodes with a run-time error in boost::math::tools::promote_args. I guess I don't really understand the purpose of promote_args. What is the purpose of promote_args? How am I using it wrong? See any other obvious problems?

Copied from original issue: stan-dev/stan/pull/1203

Update

From @bob-carpenter on September 20, 2013 20:14

mdivide_left_spd is in the code, but not exposed to Stan. Expose it and the right-handed version.

• mdivide_right_spd, defined via transpose
• manual doc
• expose function_signatures
• unit tests
• test models

Copied from original issue: stan-dev/stan#227

From @ecbrown on September 22, 2013 18:32

I am trying to implement a Matern Covariance function, which depends on K_v. I believe that v may take on non-integer values, e.g v=1/2 is quite significant from a theoretical standpoint.

However, when I supply both arguments as real (nu is a real<lower=0.5> parameter):

Sigma[i,j] <- sigmasq * 1.0/(tgamma(nu)*pow(2.0,nu-1.0) *
pow(sqrt(2.0*nu)*D[i,j]/phi,nu) * modified_bessel_second_kind(nu,
sqrt(2.0*nu)*D[i,j]/phi);

I get a compilation error:

DIAGNOSTIC(S) FROM PARSER:
no matches for function name="modified_bessel_second_kind"
   arg 0 type=real
   arg 1 type=real
available function signatures for modified_bessel_second_kind:
0.  modified_bessel_second_kind(int, real) : real

which seems to make me think that K_v only takes integer v.

Is it possible to make it take two real arguments?

Copied from original issue: stan-dev/stan#231

From @bob-carpenter on January 12, 2015 19:20

Jan Scholz (via e-mail) suggested implementing the general form of interpolation used in SciPy.

https://github.com/scipy/scipy/blob/v0.15.0/scipy/ndimage/interpolation.py#L230

It should be able to handle 2D and 3D structures. Jan says the linear form (order 1) should be enough for his purposes.

For background, BUGS provides an interp.lin function for the 1D case:

http://www.mrc-bsu.cam.ac.uk/wp-content/uploads/manual14.pdf

The goal would be to get the derivatives working properly in a way that encapsulates the index-fiddling relative to the user and takes care of error conditions.

This would be a nice standalone project.

• example for model-examples
• documentation using example
• C++ function with tests
• plumbing into Stan's language
• tests in Stan programs to ensure compilability

Copied from original issue: stan-dev/stan#1214

List of changes.

• src/stan/prob/distributions/univariate/continuous/chi_square.hpp, green 352 --- this line const T_partials_return Pn = 1.0 - gamma_p(alpha_dbl, beta_dbl * y_dbl); followed by
  ccdf_log += log(Pn); is bad. this should at least be using log1m(...).
• src/stan/prob/distributions/univariate/continuous/exp_mod_normal.hpp lots of repeated subexpressions
• same file --- this kind of line break is very confusing:

scaled_diff * SQRT_2 - deriv_2 
* sigma_dbl * SQRT_2

Code should follow the same conventions as mathematics typesetting, where this would be rendered as

scaled_diff * SQRT_2 
- deriv_2  * sigma_dbl * SQRT_2

That way, the line break reinforces the operator precedence rather than confusing it. Same thing on line green 235 (and elsewhere if you find other instances). Similarly, I would break the lines starting on green 239,

operands_and_partials.d_x4[n] += exp(0.5 * v_sq - u) 
             * (SQRT_2 / sqrt_pi * 0.5 * sigma_dbl 
                * exp(-(v / SQRT_2 - scaled_diff) * (v / SQRT_2 - scaled_diff)) 
                - (v * sigma_dbl + mu_dbl - y_dbl) * erf_calc) / cdf_;

as

operands_and_partials.d_x4[n] 
  += exp(0.5 * v_sq - u) 
     * (SQRT_2  / sqrt_pi * 0.5 * sigma_dbl 
        * exp(-square(v / SQRT_2 - scaled_diff))
        - (v * sigma_dbl + mu_dbl - y_dbl) * erf_calc)
     / cdf_;

No big deal not updating these, though I think it does make the code much easier to read.

• src/stan/prob/distributions/univariate/continuous/exp_mod_normal.hpp, green 467 --- this is another log1m situation --- maybe just searching for log(ccdf) will find them all --- it's one nice benefit of consistent naming. Oh, and I see the same thing with log(Pn) in the beta and gamma distribution code, and in 23 other places I'm not going to try to enumerate.
• src/stan/prob/distributions/univariate/continuous/exponential.hpp, green 160 and elsewhere ---The expression 1.0 - exp(...) is going to be very prone to unnecessary underflow as the inner term approaches 0, because the smallest value 1 - exp(...) can return is around 1e-14, whereas double-precision floating point should be able to accomodate closer to 1e-300. It should be more stable to return -expm1(...).
• gumbel.hpp, line 286 green --- The built-in function log1m_exp(...) is more robust than evaluating log(1 - exp(...)) directly or even log1m(exp(...)).
• green 736, grad_F32_test.cpp --- To test within 1e-5, it's not necessary to have this many digits:
  EXPECT_NEAR(0.415359887777218792995404669803015764396172842233556866773418,grad1[0],1e-5);

Originally stan-dev/stan#819 and Stan issue #706

Develop an autodiff testing framework that will test autodiff implementations against existing double implementation using signatures like:

MatrixXd x = ...;
MatrixXd y = ...;
test_ad(multiply, x, y);

The test will make sure that multiply supplies appropriate values and derivatives as defined by the double implementation and finite differences. Due to the inaccuracy of finite differences, we'll probably need an optional tolerance argument.

From @jpritikin on January 10, 2015 14:17

Per http://eigen.tuxfamily.org/dox/TopicFunctionTakingEigenTypes.html, it is more flexible if all Eigen arguments are declared as MatrixBase. This issue exists to track the conversion.

Copied from original issue: stan-dev/stan#1209

From @bob-carpenter on April 2, 2014 0:25

Kyle Foreman reports on stan-users:

...the log_determinant of the dense matrix was substantially slower in Stan than R. Last time I tried it in Stan was in December so maybe things have changed, but at that time the log determinant of a 2000x2000 matrix took on the order of 100x longer to compute in Stan than R.

I was computing it in the transformed data block, at 1000 samples of p (which is the only parameter that the log determinant varies by, and is constrained [0,1)). That way I only have to do it once (well, 1000 times, but just during initialization), and in the model block can then simply do a linear interpolation of the log determinant based on the value of parameter p. 

Maybe there's a better algorithm than what we're using? Looks like we can use Cholesky decomposition in LDLT if the matrix is positive definite.

Copied from original issue: stan-dev/stan#613

Remove build documentation from repository.

From @betanalpha on July 2, 2014 13:2

Right now the Stan autodiff implements functions with improper derivatives (i.e. the derivatives do not exist at certain points), which is formally improper. We have seen issues where people use functions like step, floor, and ceil in their models and the resulting sampling is very poor.

Personally I'd prefer to limit Stan to functions with proper derivatives only, but in any case we should have the discussion.

Copied from original issue: stan-dev/stan#731

From @bob-carpenter on October 4, 2014 16:23

Currently, stan/math.hpp directly includes these:

#include <stan/math/rep_array.hpp>
#include <stan/math/rep_vector.hpp>
#include <stan/math/rep_row_vector.hpp>
#include <stan/math/rep_matrix.hpp>

They're on the wrong path. The first (for std::vector) should go into math/functions and the remaining three (for Eigen::Matrix) into math/matrix. The tests will have to move, too.

Copied from original issue: stan-dev/stan#1067

From @bob-carpenter on February 12, 2014 18:9

It'd be nice to have the Zipf distribution in Stan:

See: http://mathworld.wolfram.com/ZipfDistribution.html

The pmf is involves the Riemann zeta function, which is available from Boost:

http://www.boost.org/doc/libs/1_55_0/libs/math/doc/html/math_toolkit/zetas/zeta.html

so we could add it directly with auto-diff. But auto-diffing these numerical
approximations to functions isn't ideal.

Plus, we really need log_zeta for Stan, and then we'd want a direct scheme
for evaluating both

log(zeta(s)

and

  d/ds log(zeta(s)) = 1/zeta(s) d/ds zeta(s).

See also, this paper by Choudhury (1995), which is what the Mathworld page cites: https://www.jstor.org/stable/52768?seq=1#page_scan_tab_contents

Copied from original issue: stan-dev/stan#551

We might want to have all of the built-in functions (from <cmath> and <math.h>) throw exceptions when they encounter domain errors, illegal arguments, NaN, underflow, or overflow.

We can't change the <cmath> and <math.h> functions themselves, but adding stan::math versions would likely cause no end of headaches for users of the autodiff library.

Originally from @mbrubake on stan-users and moved from stan-dev/stan#794

... it may be worth checking for underflow in exp() and throwing an exception. I hate to do it because it's going to potentially really hurt performance, but issues like this seem to pop up more than I'd like.

From @syclik on November 18, 2013 18:37

This really shouldn't be in the stan::agrad namespace.

Copied from original issue: stan-dev/stan#401

From @bob-carpenter on December 4, 2014 2:48

Aki Vehtari suggested building in a function that behaves like some R survival analysis packages:

real weibull_right_censored_log(real y, real sigma, real alpha, int z) {
  if (z == 0) return weibull_log(y,alpha,sigma);
  else return weibull_ccdf_log(y,alpha,sigma);
}

The int value z is used as an indicator as to whether the event was right censored at y or whether it was an ordinary outcome.

Aki also mentioned that it's common to have a vector of censored outcomes, all censored at the same place.

Copied from original issue: stan-dev/stan#1153

From @randommm on November 20, 2014 1:55

The following program:

transformed data {
  int y[2];
  y[1] <- 0;
  y[2] <- 1;
  print(poisson_ccdf_log(0, 5));
  print(poisson_ccdf_log(1, 5));
  print(poisson_ccdf_log(y, 5));
}
parameters {
  real any;
}
model {
  any ~ normal(0, 1);  
}

Outputs:

-0.00676075
-0.0412676
-0.0480283

(the last number is the sum of the first two as expected)

However, the following program:

transformed data {
  int y[2];
  y[1] <- -1;
  y[2] <- 1;
  print(poisson_ccdf_log(-1, 5));
  print(poisson_ccdf_log(1, 5));
  print(poisson_ccdf_log(y, 5));
}
parameters {
  real any;
}
model {
  any ~ normal(0, 1);  
}

Outputs:

0
-0.0412676
0

Which doesn't make sense: the last number should have been -0.0412676 too.

The way the CCDFs are currently designed makes them ignore everything when there is a negative random variable inside the vector of random variables:

  // Explicit return for extreme values
  // The gradients are technically ill-defined, but treated as neg infinity
  for (size_t i = 0; i < stan::length(n); i++) {
    if (value_of(n_vec[i]) < 0) 
      return operands_and_partials.to_var(0.0);
  }

Copied from original issue: stan-dev/stan#1142

From @betanalpha on October 11, 2014 19:7

At the moment cholesky_decompose assumes that the input matrix is positive-definite. When it is not there is no error and the output is given values as if the factorization is incorrect. Needless to say, this leads to some frustrating debugging. Either cholesky_decompose needs to throw an error if the Eigen fail bit is set, check for positive-definiteness, or require cov_matrix inputs.

Copied from original issue: stan-dev/stan#1076

From @jzelner on February 14, 2014 14:59

It would be great if the incomplete beta function was available in addition to the regular beta function.

Copied from original issue: stan-dev/stan#562

From @bob-carpenter on January 3, 2015 22:48

There's a very large (factor of around 2 at second order, more at third order) efficiency in the current implementation of fvar for our intended uses. It currently uses the template type for both values and tangents.

template <typename T>
struct fvar {
  T val_; 
  T d_;
  ...
};

This winds up propagating all kinds of derivative information to values which we never use.

New Approach

What we do use is fvar<var> and fvar<fvar<var> >. In both of these cases, we can get away with a double value type.

• Change storage type to the following
• Deal with resulting test and compilation storm

struct fvar {
  double val_; 
  T d_;
  ...
};

Example of Savings at Second Order

For example, multiplying two fvar instances a and b gives you value a.val_ * b.val_ and tangent a.d_ * b.val_ + a.val_ * b.d_.

var: value) var * var; tangent) var * var, var * var, var + var
double: value) double * double; tangent) var * double, double * var, var + var

double * double: 0 nodes, 1 multiply
var * double: 1 node, 1 multiply (val), 1 multiply (adjoint), 1 add (adjoint)
var * var: 1 node, 1 multiply (value), 2 multiplies (adjoints), 2 adds (adjoint)
var + var: 1 node, 1 add (value) 2 adds (adjoint)

Total original: 4 nodes, 9 multiply, 9 add
Total new: 3 nodes, 5 multiply, 5 add

By removing a node, there's less pointer chasing to chain() implementations and one fewer constructor called.

As an added bonus, it'll save 25% of the memory used in a multiplication (half in an addition).

Copied from original issue: stan-dev/stan#1199

From @bob-carpenter on October 9, 2014 17:13

Extend the power operator (^) to signature (matrix,int):matrix, with definitions

• A^0 = I
• A^n = A^(n-1) * A for n > 0

Throw an exception if A is not square or if n is not positive.

This should be easy to build. Extra credit if there's something clever that can be done for the gradients of each A^n[i,j] with respect to each A[i',j']. Just storing the gradients is going to be order K^4 for a K by K matrix.

Copied from original issue: stan-dev/stan#1073

From @bob-carpenter on September 23, 2014 17:19

Sebastian Weber reports that it currently allows non-positive definite covariance matrices through.

• add test in rng code and throw exception if not positive definite
• add doc in multi_normal_rng that suggests a more efficient approach based on transforming independent unit normals given a Cholesky factor

data {
  vector<lower=0>[3] sigma_eta;     // scales
  matrix[3,3] R_eta;            // correlation matrix of random effects  
}
transformed data {
  vector[3] zeros;
  zeros <- rep_vector(0.0, 3);
}
parameters {
  real<lower=0,upper=1> dummy;
}
model {
}
generated quantities {
  vector[3] eta;
  eta <- multi_normal_rng(zeros, diag_matrix(sigma_eta) * R_eta * diag_matrix(sigma_eta) );
}

and some R code to tickle it:

library(rstan)

truth <- list(
    sigma_eta = c(0.1, 0.4, 0.2),
    rho = c(0.7, 0.3, -0.5)
    )

R_eta <- diag(length(truth3$sigma_eta))/2
R_eta[1,2:3] <- truth$rho[1:2]
R_eta[2,3] <- truth$rho[3]
R_eta <- R_eta + t(R_eta)
truth$R_eta <- R_eta

## correlation matrix is not positive definite!
det(R_eta)

## but Stan does not complain at all

model <- stan("mvn-bug.stan", chains=0)

samp <- stan(fit=model, data=truth, chains=1, warmup=0, iter=1000)

p <- extract(samp, inc_warmup=FALSE, permuted=TRUE)

## now it is positive definite
Rres <- cor(p$eta)
det(Rres)

Copied from original issue: stan-dev/stan#1043

From @edbaskerville on December 9, 2014 23:33

It would be nice if Stan had some simple one-dimensional interpolation methods built in. (Not super-important, since this should be doable within a Stan model.)

One place this issue arises when implementing ODE models driven by empirical data sampled at discrete time intervals. Because the time derivative needs to be evaluated at arbitrary time points, the values of the driving variables need to be evaluated between sample times.

Even linear interpolation would be sufficient for many cases: e.g.,

data {
  int nTimepoints;
  real[nTimepoints] x;
}
// ...
x_t <- interpolate_linear(x, t);

would result in x_t equal to a value linearly interpolated from x[i], x[i + 1], where i <= t <= i + 1.

Copied from original issue: stan-dev/stan#1165

The fvar specialization for quad_form_sym doesn't currently allow for a double, fvar instantiation of the function, which makes the test/test-models/good/function-signatures/distributions/ gaussian_dlm_obs_log.hpp test fail when instantiated with fvar (or fvar for that matter)

From @bob-carpenter on December 18, 2014 21:57

• max_test.cpp: returns NaN if any input is NaN
• mdivide_left_ldlt_test.cpp
• mdivide_right_ldlt_test.cpp
• min_test.cpp
• ? fmin and fmax if they're vectorized
• sort_indices_test.cpp: if any input is NaN, then the first index will be that of the NaN (everything else that follows is what would normally be returned)
• sort_test.cpp: if any input is NaN, then the first argument in the vector return is NaN

Copied from original issue: stan-dev/stan#1176

add initial readme.

From @randommm on July 11, 2014 15:24

Add wiener first passage time density. See #469.

Copied from original issue: stan-dev/stan#765

From @bob-carpenter on November 11, 2014 16:31

• density
  • dirichlet_multinomial_lpdf function
  • distribution tests
  • model instantiation tests
  • doc for manual
  • example for manual in programming section; perhaps a new chapter on overdispersion
• PRNG
  • dirichlet_multinomial_rng function
  • distribution and PRNG tests
  • doc for manual
  • model instantiation tests

Copied from original issue: stan-dev/stan#1130

From @bob-carpenter on August 2, 2014 16:27

Joshua N. Pritikin mentioned on stan-users that the OpenMx package has a multivariate normal that works on sufficient statistics:

In OpenMx, multi_normal has two implementations. One implementation handles data in the form of 1 case per row. The other implementation handles data as a covariance matrix (and means). These two approaches are equivalent for certain data (no missingness, etc).

I would think the second implementation would be much more efficient because the sufficient statistics could be computed just once.

Does anyone know what the density over the sufficient stats would look like?

Copied from original issue: stan-dev/stan#822

Current specialization of trace_quad_form for fvar doesn't accept Eigen::Matrix<double, etc.>, Eigen::Matrix<fvar, etc.> instantiations, so multi_normal_prec_log can't be instantiated for higher-order AD.

From @syclik on November 18, 2013 18:40

src/stan/agrad/rev/matrix/Eigen_NumTraits.hpp should not be defining a class in the namespace stan::agrad. That class should be moved out of this file.

Copied from original issue: stan-dev/stan#402

From @syclik on July 8, 2014 14:25

exception handling is unlike other functions in Stan

It's probably not the only place where it's different. Should make it look like the rest of our code base.

Copied from original issue: stan-dev/stan#754

From @mbrubake on February 18, 2014 22:44

Currently still uses inverse_spd/log_determinant_spd which may be less than optimal for large systems.

Copied from original issue: stan-dev/stan#574

From @bob-carpenter on September 4, 2013 7:7

Simon Barthelmé reported via the stan-users list. I simplified the reported model to this minimal instance:

parameters {
  real x;
}
model {
  real y;
  y <- sqrt(square(x-x));

  x ~ normal(0,1);
}

As Simon correctly diagnosed, this model causes a problem with autodiff, as can be seen when running from RStan 1.3:

TESTING GRADIENT FOR MODEL 'barthelme2' NOW (CHAIN 4).

TEST GRADIENT MODE

 Log probability=-0.0571533

 param idx           value           model     finite diff           error
         0       -0.338093             nan        0.338093             nan

First step is to move this into a unit test that displays the behavior.

Copied from original issue: stan-dev/stan#203

From @bob-carpenter on November 19, 2014 19:4

To remove duplication from our derivative testing code and also make it portable for when we upgrade autodiff, we need to move our testing into a generic functional framework that takes pairs of functors and compares them.

Copied from original issue: stan-dev/stan#1141

From @bob-carpenter on November 9, 2014 6:4

It would be nice to have a way of limiting the number of steps used in the ODE integrator. The integrator could throw an exception if it hasn't converged within tolerance within a preconfigured upper bound of steps.

As is, the integrator can effectively hang while trying to achieve a desired tolerance (now hardcoded at 1e-6 for both relative and absolute).

This looks like it may be difficult with Boost odeint, at least the way we're calling it now. Maybe when we swap in an integrator that can handle stiffer systems, this problem will go away.

Copied from original issue: stan-dev/stan#1127

From @PeterLi2016 on August 8, 2014 19:39

This will speed up matrix distributions and stability of the derivatives.

For instance, all the mdivide_ functions and the trace_ functions.

Copied from original issue: stan-dev/stan#836

There are two ways to vectorize these.

Multiple I.I.D. draws

Add number of draws as first argument and generate multiple returns, e.g.,

vector[5] x = normal_rng(5, 0, 1);

Vectorized Parameters

Add vectorized versions of existing RNG functions, e.g., normal_rng(mu, sigma) where the variables mu and sigma would be either scalars or sequences such that if they were both sequences, they have the same shape. This is going to require some template metaprogramming to compute return type and also some expression templates for vector-like views of the arguments. This would allow us to have things like:

vector[3] mu = { 1, 2, 3 };
real sigma = 2.5;
vector[3] y_sim = normal_rng(mu, sigma);

A decision has to be made about return type and possible argument combinations. Would we allow normal_rng(vector, row_vector), and if so, what is the return type?

For multivariate distributions, we'd allow arrays of the input type T as std::vector<T> and we'd return std::vector<R> where R is the return type.

From @syclik on October 9, 2013 18:25

• add unit tests
• add documentation
• add unit tests for gradients

And then create an issue in stan-dev/stan to

• add to src/stan/gm/function_signatures.h
• add model instantiation tests

Copied from original issue: stan-dev/stan#271