add log_Phi_approx function and its complement

[bob-carpenter]

Description

For arithmetic stability in CDF calculations, we should have an approximate Phi and its complement on the log scale.

Note

• Ideally this will be done with custom derivatives
• Ideally this will be properly vectorized

If not, there's a log_inv_logit that could be called.

A naming convention for complements needs to be established. logc_Phi_approx seems a bit too cryptic despite our having lcdf and lccdf functions.

Current Math Version

v2.18.0

[bbbales2]

Silly trick, but could we make a few scalar subclasses and automate a lot of this log scale stuff? It wouldn't fix everything, but it could get some I think.

Like, call it sum_exp_var, and we basically overload operator+ and division. And any time we exp a var it turns into a sum_exp_var.

template<typename ...Targs>
sum_exp_var {
  std::tuple<Targs...> t;

  sum_exp_var operator+=(const sum_exp_var& other) {
    return sum_exp_var(std::tuple_cat(this->t, other.t)); // Wait till we have all the terms so can extract common factors
    // Or if we end up with a tuple == (1.0, var) or whatever we could treat that specially?
  }

  sum_exp_var operator/=(const sum_exp_var& other) {
    // Can pull out common denominators here
  }
}

And make a var constructor from sum_exp_vars, so it can happen implicitly.

And there could be specializations for log1p, log1m, and log and such?

I know there's some code for the reverse mode specializations that looks like:

enable_if<std::is_same<T, stan::math::var>, stan::math::var>::type
func(const T& t) {
  ...
}

And I think the exp_vars might screw up the intended action there (and the exp_vars end up using the prim implementations of whatever functions). Bleh. Honestly the sounds awful to debug, but it's super annoying trying to write the numerically correct versions of these things.

[bbbales2]

return sum_exp_var(std::tuple_cat(this->t, other.t));

Bah this feels like it breaks the whole knowing tuple sizes at runtime thing. Presumably it wouldn't actually compile. And if the memory wasn't fancy on these things it might end up just being super slow.

[bob-carpenter]

Yes, I think we could use expression templates to express things on the linear scale and then have them auto-transformed to the log scale in some cases. In others, it'd require some fancy algebra.

The probabilistic programming system Hakaru expresses everything on the linear scale then computes on the log scale.