sd() problem: sqrt(square(0))

[Daniel Lee]

We're getting nan's in sd() when all the elements are the same.

It's really coming from taking the sqrt of the square of 0. A minimal test case:

#include <stan/agrad/rev/sqrt.hpp>
#include <gtest/gtest.h>
#include <stan/agrad/var.hpp>

TEST(AgradRev,sqrt_square_zero) {
  stan::agrad::var x(0.0);
  stan::agrad::var f = sqrt(square(x));
 
  EXPECT_FLOAT_EQ(0.0,f.val());
  std::vector<stan::agrad::var> g;
  g.push_back(x);
  std::vector<double> gradient;
  f.grad(g,gradient);
  EXPECT_FLOAT_EQ(1.0,gradient[0]);
}

The gradient ends up being nan. By the way, this happens with multiplying x by itself as well. And if we change the function order.

The two chain methods (just the relevant pieces pulled out):

sqrt:
        void chain() {
          avi_->adj_ += adj_ / (2.0 * val_);
        }

square:
        void chain() {
          avi_->adj_ += adj_ * 2.0 * avi_->val_;
        }

In the test case outlined above, we start with the chain method in sqrt with:

avi_->adj_ = 0

adj_ = 1

avi_->val_ = inf

The resulting adjoint value is as expected: inf. (this seems reasonable in the limit as x approaches 0.)

Then we apply the chain rule and get into the chain method in sqrt with:

avi_->adj_ = 0

adj_ = inf

avi_->val_ = 0

The resulting adjoint value blows up: nan. This is from the inf * 2 * 0 condition. inf * 0 = nan.

The obvious right answer is d/dx sqrt(square(x)) = d/dx x = 1, but I can't figure out how to get there from local chain functions without some knowledge of the expression tree. An inf or 0 could arise in lots of ways, so it doesn't seem like creating a special case is the right thing to do.

So... what should we do?

Daniel

[Bob Carpenter]

I can fix it. I worked out the analytic gradients,
which don't have this problem. That'll solve the
specific case of sd().

Not much we can do about general code that does sqrt(0)
but has some analytic reduction, like sqrt(square(x)) = x.
In general, there are lots of ways to express square(x), such
as x * x. And the case for sd() is more like sqrt(sum(square(x)).

I don't know how to do runtime type resolution, but we
might be able to catch the expression tree structure and
collapse it. But that doesn't solve the multiple ways
to express it problem.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.

[Daniel Lee]

On Tue, Nov 19, 2013 at 4:50 PM, Bob Carpenter <ca...@alias-i.com> wrote:

I can fix it.  I worked out the analytic gradients,
which don't have this problem.  That'll solve the
specific case of sd().

Oh. The analytic gradients for sd() itself? That makes sense. Cool.

 

Not much we can do about general code that does sqrt(0)
but has some analytic reduction, like sqrt(square(x)) = x.
In general, there are lots of ways to express square(x), such
as x * x.  And the  case for sd() is more like sqrt(sum(square(x)).

Yes. I was trying to reduce it to the simplest case.

 

I don't know how to do runtime type resolution, but we
might be able to catch the expression tree structure and
collapse it.   But that doesn't solve the multiple ways
to express it problem.

Hmmm.

 

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

For more options, visit https://groups.google.com/groups/opt_out.

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

[Bob Carpenter]

On 11/19/13, 4:57 PM, Daniel Lee wrote:

> On Tue, Nov 19, 2013 at 4:50 PM, Bob Carpenter <ca...@alias-i.com <mailto:ca...@alias-i.com>> wrote:
>
> I can fix it. I worked out the analytic gradients,
> which don't have this problem. That'll solve the
> specific case of sd().
>
>
> Oh. The analytic gradients for sd() itself? That makes sense. Cool.

Yes, the gradient of var(x) for a vector x is a lot of chopping,
but it reduces to a very simple form. Hence, so does sd(x) and
given that mean(x) is also trivial, so does standardize(x).

I'll add the code for fwd and rev to do it analytically and
add some tests for the zero case.

- Bob

[Marcus Brubaker]

Didn't this issue come up a while back?

Mathematically, the derivative d/dx (sqrt(x)) is not well defined at x=0 so in some sense we can put whatever value we want because basically everything is wrong.  From a functional perspective, we never want to get a non-finite gradient out of the autodiff unless we absolutely must (e.g., someone divides by zero).  Putting those two things together, maybe we can just rewrite chain for sqrt() to be:

        void chain() {

          if (val_ > 1e-300)

              avi_->adj_ += adj_ / (2.0 * val_);

          else

              avi_->adj_ += adj_ / 2.0;

        }

or, alternatively,

        void chain() {

          avi_->adj_ += adj_ / (2.0 * std::max(val_,2e-300));

        }

Cheers,

Marcus

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

[Bob Carpenter]

On 11/19/13, 5:12 PM, Marcus Brubaker wrote:
> Didn't this issue come up a while back?

Yes, but we never resolved it.

This time it's a bit different because we really
only need to differentiate sd(x).

> Mathematically, the derivative d/dx (sqrt(x)) is not well defined at x=0 so in some sense we can put whatever value we
> want because basically everything is wrong. From a functional perspective, we never want to get a non-finite gradient
> out of the autodiff unless we absolutely must (e.g., someone divides by zero).

Taking a derivative of sqrt(x) where x = 0 is essentially
dividing by 0.

It hasn't come up as an issue except where someone literally
wrote sqrt(square(x)) in code as an example and in sd().

> Putting those two things together, maybe
> we can just rewrite chain for sqrt() to be:
>
> void chain() {
> if (val_ > 1e-300)
> avi_->adj_ += adj_ / (2.0 * val_);
> else
> avi_->adj_ += adj_ / 2.0;
> }
>
> or, alternatively,
>
> void chain() {
> avi_->adj_ += adj_ / (2.0 * std::max(val_,2e-300));
> }

Will this do the right thing for HMC and BFGS?

I'm also OK with just raising an exception or returning an infinite
gradient.

- Bob

[Marcus Brubaker]

This is effectively returning a subgradient for points where the gradient isn't defined.  HMC won't care so long as same subgradient is always returned (which it is in this case).  The regions where this will happen are extremely small, so, except for a few pathological cases, I wouldn't expect it to impact either method significantly.

I'm also OK with just raising an exception or returning an infinite
gradient.

That could also work, but in some sense it just moves the problem around a bit in terms of things like initialization.

Cheers,

Marcus

[Bob Carpenter]

OK -- I'll make the fix you (Marcus) suggest for
sqrt() in a separate pull request from the
sd() fix. I'll associate the sqrt() fix with
issue 203, which is 3 months old:

https://github.com/stan-dev/stan/issues/203

- Bob

On 11/19/13, 8:19 PM, Marcus Brubaker wrote:

[Michael Betancourt]

Sorry for the late reply - usual conference business.

I disagree with using this fix as HMC doesn't make sense when the potential isn't smooth (recall that sqrt(square(x)) is just the absolute value which has a cusp as zero). As has already been noted the probability of hitting the cusp is very small, but it nonzero and at a point where users often initialize. It IS an initialization problem.

I would much prefer to see an exception thrown so that users are aware of the issue. Hiding it by saying "we can define it to be whatever we want" (which is false) is bad statistical mojo,

[Bob Carpenter]

OK, is it OK with everyone to throw an exception?
Basically, returning NaN is tantamount to doing just
that in the context of Stan's log prob and gradient
evals.

We do the same thing with abs(x) itself, taking its
derivative to be 0 at 0. This seems more sensible than for
sqrt, as the limit from the left is -1 and from the right is 1.

There's a big problem with inits of 0, which can
inadvertently land you right on one of these discontinuities.

- Bob

[Michael Betancourt]

That works for me. It really should be a proper exception that let's the user know there's a problem with their model rather than sweeping it under the rug with a sampler rejection

To be fair, only the absolute value has a cusp within its support. Just about everything else has an undefined gradient at a boundary which gets mapped to infinity where we should never reach anyways.

[Ben Goodrich]

BTW, you are going to have the same issue whenever we reimplement vectors that are uniform on a sphere by normalizing a vector of iid standard normals. If they do init = 0, then you divide by the square root of zero in the denominator.

Ben

[Marcus Brubaker]

Returning a subgradient won't affect HMC at all except maybe on the level of efficiency but then only when there is significant mass at that region and even then I expect the impact to be small because we're returning something which will be close.  It's not "bad statistical mojo" as It certainly won't impact the statistical validity of our samples since a) the Metropolis test will use the same "correct" log prob and b) so long as the same subgradient is returned, the simulation will be reversible.  Patching the gradient is basically just using a different simulation Hamiltonian which is only different from the correct one in a small region around these gradient singularities.

If we throw an exception we're just going to reject everything and stop sampling in that region altogether.   Further, these rejections won't just occur at exactly x = 0, returning an unusually large gradient can overflow chain() calculations elsewhere resulting in a nan or inf gradient.

Regarding efficiency of patching the gradient vs throwing an exception, it's not clear which will be more efficient.  It seems to me allowing it to traverse through a region instead of automatically rejecting any trajectory which comes near the singularity would be more statistically efficient.

[Marcus Brubaker]

On Wed, Nov 20, 2013 at 11:02 AM, Ben Goodrich <goodri...@gmail.com> wrote:

BTW, you are going to have the same issue whenever we reimplement vectors that are uniform on a sphere by normalizing a vector of iid standard normals. If they do init = 0, then you divide by the square root of zero in the denominator.

Yeah, thankfully that should only be a problem at initialization.  It's unlikely to ever be seen during simulation.

On that tangent, I've been doing some work looking at parameterizing the sphere for sampling and it looks like there may be better options than the vector of iid normals which has zero probability density at the origin, so the problem should never arise.

Cheers,

Marcus

[Bob Carpenter]

The double-exponential (Laplace) distro uses abs(), so it
also has a cusp. It's a reasonble prior for regression
coeffs and 0 is a reasonable init. So we want this to work
and not require a slightly jittered init.

- Bob

[Daniel Lee]

On Wed, Nov 20, 2013 at 11:52 AM, Marcus Brubaker <marcus.a...@gmail.com> wrote:

Returning a subgradient won't affect HMC at all except maybe on the level of efficiency but then only when there is significant mass at that region and even then I expect the impact to be small because we're returning something which will be close.  

Unfortunately, we're artificially giving mass in these regions significant mass since we're initializing starting at 0.

 

It's not "bad statistical mojo" as It certainly won't impact the statistical validity of our samples since a) the Metropolis test will use the same "correct" log prob and b) so long as the same subgradient is returned, the simulation will be reversible.  Patching the gradient is basically just using a different simulation Hamiltonian which is only different from the correct one in a small region around these gradient singularities.

If we throw an exception we're just going to reject everything and stop sampling in that region altogether.   Further, these rejections won't just occur at exactly x = 0, returning an unusually large gradient can overflow chain() calculations elsewhere resulting in a nan or inf gradient.

Regarding efficiency of patching the gradient vs throwing an exception, it's not clear which will be more efficient.  It seems to me allowing it to traverse through a region instead of automatically rejecting any trajectory which comes near the singularity would be more statistically efficient.

I'm in favor of catching errors early. One reason this was difficult to track down from the outside was due to the gradients going to nan without informative error messages. It was baffling because it only happened with sampling and not with optimization or test gradient mode. I believe we'll get a handle on the sqrt(square(0)) case soon, but the fundamental reason for this issue isn't going away and I'd like to be able to identify the problems much quicker.

Right now, I don't think we check gradients for nan values. Should we? It seems like a reasonable thing to do as opposed to what we currently do now, which is to integrate using nan derivatives and failing when the nans propagate to the positions. This could be done independent of error checking using exceptions.

[Marcus Brubaker]

On Wed, Nov 20, 2013 at 12:12 PM, Daniel Lee <bea...@alum.mit.edu> wrote:

On Wed, Nov 20, 2013 at 11:52 AM, Marcus Brubaker <marcus.a...@gmail.com> wrote:

Returning a subgradient won't affect HMC at all except maybe on the level of efficiency but then only when there is significant mass at that region and even then I expect the impact to be small because we're returning something which will be close.  

Unfortunately, we're artificially giving mass in these regions significant mass since we're initializing starting at 0.

That doesn't give the region mass.  If it has little real mass, sampling should quickly transition out of the region assuming that sampling can actually proceed.

 

 

It's not "bad statistical mojo" as It certainly won't impact the statistical validity of our samples since a) the Metropolis test will use the same "correct" log prob and b) so long as the same subgradient is returned, the simulation will be reversible.  Patching the gradient is basically just using a different simulation Hamiltonian which is only different from the correct one in a small region around these gradient singularities.

If we throw an exception we're just going to reject everything and stop sampling in that region altogether.   Further, these rejections won't just occur at exactly x = 0, returning an unusually large gradient can overflow chain() calculations elsewhere resulting in a nan or inf gradient.

Regarding efficiency of patching the gradient vs throwing an exception, it's not clear which will be more efficient.  It seems to me allowing it to traverse through a region instead of automatically rejecting any trajectory which comes near the singularity would be more statistically efficient.

I'm in favor of catching errors early. One reason this was difficult to track down from the outside was due to the gradients going to nan without informative error messages. It was baffling because it only happened with sampling and not with optimization or test gradient mode. I believe we'll get a handle on the sqrt(square(0)) case soon, but the fundamental reason for this issue isn't going away and I'd like to be able to identify the problems much quicker.

Catching errors early is good but throwing exceptions is only useful if the caller is actually in a position to do something useful about it or there is no other resory.  In this case, all the caller can do is reject the sample and/or bail back to the user who is even less equipped to do something about it.

 

Right now, I don't think we check gradients for nan values. Should we? It seems like a reasonable thing to do as opposed to what we currently do now, which is to integrate using nan derivatives and failing when the nans propagate to the positions. This could be done independent of error checking using exceptions.

We should definitely be checking gradients for non-finite values.

Cheers,

Marcus

[Michael Betancourt]

I thought a NaN gradient would throw the same exception that bad constraints does
but I guess I'm wrong. We could just have the gradient functors check and throw
the necessary exceptions.

[Daniel Lee]

On Wed, Nov 20, 2013 at 12:24 PM, Marcus Brubaker <marcus.a...@gmail.com> wrote:

On Wed, Nov 20, 2013 at 12:12 PM, Daniel Lee <bea...@alum.mit.edu> wrote:

On Wed, Nov 20, 2013 at 11:52 AM, Marcus Brubaker <marcus.a...@gmail.com> wrote:

Returning a subgradient won't affect HMC at all except maybe on the level of efficiency but then only when there is significant mass at that region and even then I expect the impact to be small because we're returning something which will be close.  

Unfortunately, we're artificially giving mass in these regions significant mass since we're initializing starting at 0.

That doesn't give the region mass.  If it has little real mass, sampling should quickly transition out of the region assuming that sampling can actually proceed.

Our problem lies with implementation. I agree with everything you've said in theory, but the implementation has initialization start in that region. We can't and don't get out. Stan just comes to a halt with an error message stating:

 

 

 

It's not "bad statistical mojo" as It certainly won't impact the statistical validity of our samples since a) the Metropolis test will use the same "correct" log prob and b) so long as the same subgradient is returned, the simulation will be reversible.  Patching the gradient is basically just using a different simulation Hamiltonian which is only different from the correct one in a small region around these gradient singularities.

If we throw an exception we're just going to reject everything and stop sampling in that region altogether.   Further, these rejections won't just occur at exactly x = 0, returning an unusually large gradient can overflow chain() calculations elsewhere resulting in a nan or inf gradient.

Regarding efficiency of patching the gradient vs throwing an exception, it's not clear which will be more efficient.  It seems to me allowing it to traverse through a region instead of automatically rejecting any trajectory which comes near the singularity would be more statistically efficient.

I'm in favor of catching errors early. One reason this was difficult to track down from the outside was due to the gradients going to nan without informative error messages. It was baffling because it only happened with sampling and not with optimization or test gradient mode. I believe we'll get a handle on the sqrt(square(0)) case soon, but the fundamental reason for this issue isn't going away and I'd like to be able to identify the problems much quicker.

Catching errors early is good but throwing exceptions is only useful if the caller is actually in a position to do something useful about it or there is no other resory.  In this case, all the caller can do is reject the sample and/or bail back to the user who is even less equipped to do something about it.

I'm assuming the error messages would bubble up past the caller all the way up to the user. Or at least to a developer that's trying to find issues with bad gradients.

[Daniel Lee]

Sorry, the current error message is:
"Rejecting inititialization at zero because of divergent gradient."

Here's the model (sd.stan):
parameters {
  vector<lower=-1,upper=1>[5] y;
}
transformed parameters {
  real sd_y;
  sd_y <- sd(y);
}
model {
  y ~ normal(0,1);
}

Here's the command line output:
> ../user_models/issue-379/sd sample init=0
Rejecting inititialization at zero because of divergent gradient.

>

That's our current implementation. We start there and end there. No indication that the sd() function is the one that's giving us grief.

[Marcus Brubaker]

On Wed, Nov 20, 2013 at 12:45 PM, Daniel Lee <bea...@alum.mit.edu> wrote:

Sorry, the current error message is:
"Rejecting inititialization at zero because of divergent gradient."

Here's the model (sd.stan):
parameters {
  vector<lower=-1,upper=1>[5] y;
}
transformed parameters {
  real sd_y;
  sd_y <- sd(y);
}
model {
  y ~ normal(0,1);
}

Here's the command line output:
> ../user_models/issue-379/sd sample init=0
Rejecting inititialization at zero because of divergent gradient.

>

That's our current implementation. We start there and end there. No indication that the sd() function is the one that's giving us grief.

If we change the gradient near singularities does it still get stuck?  Throwing an exception effectively punts the problem to the user.  I agree that we should be producing better error better error messages.

[Daniel Lee]

On Wed, Nov 20, 2013 at 12:38 PM, Michael Betancourt <betan...@gmail.com> wrote:

I thought a NaN gradient would throw the same exception that bad constraints does
but I guess I'm wrong.  We could just have the gradient functors check and throw
the necessary exceptions.

No, the checks on probability distributions only check the values, not adjoints. And the basic math functions don't have the same checks. We behave like the standard template library and default to their behavior on taking in nans, infs, 0s, and other boundary conditions where reasonable.

 

[Daniel Lee]

On Wed, Nov 20, 2013 at 12:48 PM, Marcus Brubaker <marcus.a...@gmail.com> wrote:

If we change the gradient near singularities does it still get stuck?

No, it moves and samples just fine.

 

 Throwing an exception effectively punts the problem to the user.  I agree that we should be producing better error better error messages.

Yes, exactly. It's going to be difficult to hunt down every singular point. If users are alerted to a problem, if it's an actual problem with Stan code, I'm trusting that enough people will bring it to our attention. So yes, punting the problem to the user seems like the right thing to do.

[Michael Betancourt]

There are two options -- either there's not enough mass around zero, chains never hit the singularity, and everything is fine. But if there is enough mass around zero then the chains can hit the singularity and problems arise. As Dan noted the fact that people like "0" makes initialization a problem even if sampling isn't. In either case, the user should be informed because they have an ill-defined gradient being used in their model which is a bad idea™ as it is usually indicative of a modeling (or initialization) mistake.

[Daniel Lee]

On Wed, Nov 20, 2013 at 12:56 PM, Michael Betancourt <betan...@gmail.com> wrote:

There are two options -- either there's not enough mass around zero, chains never hit the singularity, and everything is fine.  But if there is enough mass around zero then the chains can hit the singularity and problems arise.  As Dan noted the fact that people like "0" makes initialization a problem even if sampling isn't.  In either case, the user should be informed because they have an ill-defined gradient being used in their model which is a bad idea™ as it is usually indicative of a modeling (or initialization) mistake.

Speaking of initialization mistakes.... is init_adapt() supposed to start with all the unconstrained parameters at 0? (I don't recall if this ever got answered)

Here's the relevant code from src/stan/gm/command.hpp:
            case 111: {
              typedef stan::mcmc::adapt_diag_e_nuts<Model, rng_t> sampler;
              sampler_ptr = new sampler(model, base_rng, num_warmup, &std::cout, &std::cout);
              if (!init_nuts<sampler>(sampler_ptr, algo)) return 0;
              if (!init_adapt<sampler>(sampler_ptr, adapt)) return 0;
              break;
            }

Notice, the initialization, the cont_params, never gets used here. So it always starts at 0.

 

On Nov 20, 2013, at 6:48 PM, Marcus Brubaker <marcus.a...@gmail.com> wrote:

> If we change the gradient near singularities does it still get stuck?  Throwing an exception effectively punts the problem to the user.  I agree that we should be producing better error better error messages.

[Michael Betancourt]

init_adapt is just setting the adaption parameters (either to their defaults or to the user-specified values) -- it doesn't touch the initial state of the sampler.

[Daniel Lee]

Unfortunately, no. Inside init_adapt(), line 231, "s->init_stepsize();" has some side-effects and this is where the current code is bombing.

This is calling the init_stepsize() function in src/stan/mcmc/base_hmc.hpp and this is where things are initialized to 0.

[Marcus Brubaker]

But what's our response once it's been brought to our attention?  Should we work around it by patching gradients where reasonable?

Anyway, just to put a clear plan out there for us to discuss:

1) We should be checking gradients for non-finite values and reporting this as an error as soon as it's reasonable to detect.

2) Reporting more detailed error messages would be nice if we can find a way to do it without killing performance.  Perhaps we could have a compile-time option to check all adjoints after every individual chain() call.  This will be expensive, but will tell us precisely where our gradients go awry.  Alternatively we could add code to every chain() function but that's probably more work than simply hunting down every singular point.

3) When such errors are found and reported I think we should patch the gradients where it's reasonable to do so.  If the function itself is well defined at a point, we can generally come up with a reasonable subgradient which will allow sampling to proceed instead of simply punting the problem up the chain where there is less ability to do anything productive about it.

Basically all of the unconstrained parameters have been defined so that 0 inits should provide a reasonable starting point.  If a reasonable model (e.g., the double exponential) cannot be initialized at zero that seems like something we should fix as opposed to just saying it's a modelling mistake or requiring our users to jump through the massive hoops required for an explicit manual initialization.  I think this would be a somewhat different discussion if we had manual initialization block so we could just tell the user "change your initialization for this variable".

Cheers,

Marcus

[Bob Carpenter]

I think this boils down to what we'd like to
see happen for the following calculation where
x and y are parameters:

mean <- (x + y) / 2;
var <- (pow(x-mean,2) + pow(y-mean,2)) / 2;
sd <- sqrt(var);

There's nothing wrong with the mathematical definition
and the variable sd has a finite analytic gradient w.r.t.
x and y.

The problem is that the current auto-diff's brute-force application
of the chain rule is going to kill the code on the sqrt(0)
if x == y.

We can write our own sd() with analytic gradients and solve
the problem for that case.

But the real issue is what to do for the above model?

As a side note, we should be promoting diagnose mode (test_grad
in RStan more heavily.

- Bob

On 11/20/13, 12:56 PM, Michael Betancourt wrote:
> There are two options -- either there's not enough mass around zero, chains never hit the singularity, and everything is fine. But if there is enough mass around zero then the chains can hit the singularity and problems arise. As Dan noted the fact that people like "0" makes initialization a problem even if sampling isn't. In either case, the user should be informed because they have an ill-defined gradient being used in their model which is a bad idea� as it is usually indicative of a modeling (or initialization) mistake.

[Andrew Gelman]

Hi, I'm confused. Var should divide by n-1, not by n.
I'm not joking here, dividing by n-1 is the definition.
We can have the convention that if n=1, we define variance to be 0.

On Nov 20, 2013, at 7:48 PM, Bob Carpenter wrote:

> I think this boils down to what we'd like to
> see happen for the following calculation where
> x and y are parameters:
>
> mean <- (x + y) / 2;
> var <- (pow(x-mean,2) + pow(y-mean,2)) / 2;
> sd <- sqrt(var);
>
> There's nothing wrong with the mathematical definition
> and the variable sd has a finite analytic gradient w.r.t.
> x and y.
>
> The problem is that the current auto-diff's brute-force application
> of the chain rule is going to kill the code on the sqrt(0)
> if x == y.
>
> We can write our own sd() with analytic gradients and solve
> the problem for that case.
>
> But the real issue is what to do for the above model?
>
> As a side note, we should be promoting diagnose mode (test_grad
> in RStan more heavily.
>
> - Bob
>
> On 11/20/13, 12:56 PM, Michael Betancourt wrote:

>> There are two options -- either there's not enough mass around zero, chains never hit the singularity, and everything is fine. But if there is enough mass around zero then the chains can hit the singularity and problems arise. As Dan noted the fact that people like "0" makes initialization a problem even if sampling isn't. In either case, the user should be informed because they have an ill-defined gradient being used in their model which is a bad idea™ as it is usually indicative of a modeling (or initialization) mistake.

>>
>> On Nov 20, 2013, at 6:48 PM, Marcus Brubaker <marcus.a...@gmail.com> wrote:
>>
>>> If we change the gradient near singularities does it still get stuck? Throwing an exception effectively punts the problem to the user. I agree that we should be producing better error better error messages.
>>
>

[Bob Carpenter]

[pulled into its own thread for semantics mavens]

On 11/20/13, 1:50 PM, Andrew Gelman wrote:
> Hi, I'm confused. Var should divide by n-1, not by n.
> I'm not joking here, dividing by n-1 is the definition.

Sorry --- bad example. Stan absolutely divides by N-1
in the sd() and var() functions.

I'd go even further and separate four concepts that are
often just called "variance" in the literature:

1. variance (of a random variable or distribution)

2. population variance

3. sample variance

4. estimate of variance based on a sample

For (1), everyone divides by N for finite discrete random variables

For (2), Snedecor and Cochran (closest reliable book at hand) define
population variance as dividing by N.

For (3), Snedecor and Cochran define sample variance as dividing by N-1.

For (4), we know there are many estimators, the MLE (divide by N)
and unbiased (divide by N-1) being the usual pair in intro math
stats books. Arguably, a point estimate of variance should not
be called "variance", no matter how you define it.

Yet one more reason intro stat students (like I was only very very
recently and arguably still am) get really
really really confused from the standard presentations of this
material in intro courses and books and tutorials. For instance,
the R doc, always a source of enlightment, says simply:

This function computes the standard deviation of the values in x.

and later in the details section:

Like var this uses denominator n - 1.

The standard deviation of a zero-length vector (after removal of
NAs if na.rm = TRUE) is not defined and gives an error. The standard deviation
of a length-one vector is NA.

> We can have the convention that if n=1, we define variance to be 0.

We return the special NaN value rather than returning 0. We
treat it as 0 in n_eff calcs, though, I believe, so we don't get
all those NaN cases when a param doesn't move in a chain.

- Bob

[Andrew Gelman]

I'd set sd of a length 1 vector to the value 0. No point in having Nans running around, and I can easily envision examples where I might want to take sd in general cases that might include length 1 vectors.
The way it would work is that we might compute a weighted avg with a term such as (n-1)*sd. So the n-1 zeroes it out anyway. I don't think we should be giving out Nans when we can avoid it.

[Bob Carpenter]

Sorry --- I keep confusing the 0 and 1 cases.
If size(x) == 0, then we raise an error. We could
have this case also return 0 and thus deal with the boundary
condition of empty bins in models.

In Stan 2.0.1, if size(x) == 1, then

sd(x) = var(x) = 0

and

d/dx_1 sd(x) = d/dx_1 var(x) = 0

Whatever we decide, I'll improve the doc in the manual to reflect
the boundary cases. We really need to do this for all of our functions,
as it's coming up a lot in user models. So expect about five
dozen of these discussions for everything from trig to Bessel
functions :-)

- Bob