Make normal_lcdf(x | 0, 1) equivalent to Phi()

@PhilClemson

Since #1411 normal_lcdf is more robust in the tails than Phi or even Phi_approx. Update Phi() to be as robust as normal_lcdf and add log_Phi() as discussed in the linked PR above (should this be a separate issue?)

In fact, Phi_approx() is less robust than Phi(). Adding @PhilClemson @nhuurre @bbbales2

generated quantities {
  vector[N] b;
  vector[N] c;
  vector[N] d;
   for (i in 1:N) {
    b[i] = normal_lcdf(i | 0, 1);
    c[i] = log(Phi_approx(i));
    d[i] = log(Phi(i));
   }
}

> data.table(mean_b = test_b, mean_c = test_c, mean_d = test_d)
           mean_b       mean_c       mean_d
 1:  -1.72754e-01 -1.72771e-01 -1.72754e-01
 2:  -2.30129e-02 -2.30241e-02 -2.30129e-02
 3:  -1.35081e-03 -1.23271e-03 -1.35081e-03
 4:  -3.16717e-05 -1.83432e-05 -3.16717e-05
 5:  -2.86652e-07 -5.01624e-08 -2.86652e-07
 6:  -9.86588e-10 -1.65181e-11 -9.86588e-10
 7:  -1.27981e-12 -4.44089e-16 -1.27987e-12
 8:  -6.22096e-16  0.00000e+00 -6.66134e-16
 9:  -1.12859e-19  0.00000e+00  0.00000e+00
10:  -7.61985e-24  0.00000e+00  0.00000e+00
11:  -1.91066e-28  0.00000e+00  0.00000e+00
12:  -1.77648e-33  0.00000e+00  0.00000e+00
13:  -6.11716e-39  0.00000e+00  0.00000e+00
14:  -7.79354e-45  0.00000e+00  0.00000e+00
15:  -3.67097e-51  0.00000e+00  0.00000e+00
16:  -6.38875e-58  0.00000e+00  0.00000e+00
17:  -4.10600e-65  0.00000e+00  0.00000e+00
18:  -9.74095e-73  0.00000e+00  0.00000e+00
19:  -8.52722e-81  0.00000e+00  0.00000e+00
20:  -2.75362e-89  0.00000e+00  0.00000e+00
21:  -3.27928e-98  0.00000e+00  0.00000e+00
22: -1.43989e-107  0.00000e+00  0.00000e+00
23: -2.33064e-117  0.00000e+00  0.00000e+00
24: -1.39039e-127  0.00000e+00  0.00000e+00
25: -3.05670e-138  0.00000e+00  0.00000e+00
26: -2.47606e-149  0.00000e+00  0.00000e+00
27: -7.38948e-161  0.00000e+00  0.00000e+00
28: -8.12387e-173  0.00000e+00  0.00000e+00
29: -3.28979e-185  0.00000e+00  0.00000e+00
30: -4.90671e-198  0.00000e+00  0.00000e+00
31: -2.69525e-211  0.00000e+00  0.00000e+00
32: -5.45208e-225  0.00000e+00  0.00000e+00
33: -4.06119e-239  0.00000e+00  0.00000e+00
34: -1.11390e-253  0.00000e+00  0.00000e+00
35: -1.12491e-268  0.00000e+00  0.00000e+00
36: -4.18262e-284  0.00000e+00  0.00000e+00
37: -5.72557e-300  0.00000e+00  0.00000e+00
38: -2.88543e-316  0.00000e+00  0.00000e+00
39:   0.00000e+00  0.00000e+00  0.00000e+00
40:   0.00000e+00  0.00000e+00  0.00000e+00

It looks while log(Phi(x)) does not perform as well as std_normal_lcdf(x), there is only a marginal difference between Phi(x) and exp(std_normal_lcdf(x)). Also, log(Phi_approx(x)) performs better (but not equivalently) when log_inv_logit is used.

Testing

Model functions

funs = "
functions {
  real log_phi_old(real x) {
    return log(Phi(x));
  }
  real log_phi_new(real x) {
    return std_normal_lcdf(x|);
  }
  real log_phi_approx_old(real x) {
    return log(Phi_approx(x));
  }
  real log_phi_approx_new(real x) {
    return log_inv_logit(0.07056 * pow(x, 3.0) + 1.5976 * x);
  }
  real phi_old(real x) {
    return Phi(x);
  }
  real phi_new(real x) {
    return exp(std_normal_lcdf(x|));
  }
  real phi_approx_old(real x) {
    return Phi_approx(x);
  }
  real phi_approx_new(real x) {
    return exp(log_inv_logit(0.07056 * pow(x, 3.0) + 1.5976 * x));
  }
}
"

rstan::expose_stan_functions(rstan::stanc(model_code=funs))

`log(Phi(x))` Results:

get_logphi = function(x) {
    data.frame(
        x = x,
        LPhi_new = log_phi_new(x),
        LPhi_old = log_phi_old(x),
        LPhi_appr_new = log_phi_approx_new(x),
        LPhi_appr_old = log_phi_approx_old(x)
    )
}

purrr::map_dfr(1:40,get_logphi)

    x       LPhi_new      LPhi_old  LPhi_appr_new LPhi_appr_old
1   1  -1.727538e-01 -1.727538e-01  -1.727709e-01 -1.727709e-01
2   2  -2.301291e-02 -2.301291e-02  -2.302409e-02 -2.302409e-02
3   3  -1.350810e-03 -1.350810e-03  -1.232715e-03 -1.232715e-03
4   4  -3.167174e-05 -3.167174e-05  -1.834324e-05 -1.834324e-05
5   5  -2.866516e-07 -2.866516e-07  -5.016240e-08 -5.016240e-08
6   6  -9.865876e-10 -9.865877e-10  -1.651817e-11 -1.651812e-11
7   7  -1.279813e-12 -1.279865e-12  -4.289120e-16 -4.440892e-16
8   8  -6.220961e-16 -6.661338e-16  -5.750875e-22  0.000000e+00
9   9  -1.128588e-19  0.000000e+00  -2.607333e-29  0.000000e+00
10 10  -7.619853e-24  0.000000e+00  -2.617536e-38  0.000000e+00
11 11  -1.910660e-28  0.000000e+00  -3.810306e-49  0.000000e+00
12 12  -1.776482e-33  0.000000e+00  -5.266657e-62  0.000000e+00
13 13  -6.117164e-39  0.000000e+00  -4.526424e-77  0.000000e+00
14 14  -7.793537e-45  0.000000e+00  -1.584009e-94  0.000000e+00
15 15  -3.670966e-51  0.000000e+00 -1.478016e-114  0.000000e+00
16 16  -6.388754e-58  0.000000e+00 -2.408003e-137  0.000000e+00
17 17  -4.105996e-65  0.000000e+00 -4.485680e-163  0.000000e+00
18 18  -9.740949e-73  0.000000e+00 -6.256481e-192  0.000000e+00
19 19  -8.527224e-81  0.000000e+00 -4.278579e-224  0.000000e+00
20 20  -2.753624e-89  0.000000e+00 -9.394498e-260  0.000000e+00
21 21  -3.279278e-98  0.000000e+00 -4.337000e-299  0.000000e+00
22 22 -1.439892e-107  0.000000e+00   0.000000e+00  0.000000e+00
23 23 -2.330637e-117  0.000000e+00   0.000000e+00  0.000000e+00
24 24 -1.390392e-127  0.000000e+00   0.000000e+00  0.000000e+00
25 25 -3.056697e-138  0.000000e+00   0.000000e+00  0.000000e+00
26 26 -2.476063e-149  0.000000e+00   0.000000e+00  0.000000e+00
27 27 -7.389481e-161  0.000000e+00   0.000000e+00  0.000000e+00
28 28 -8.123869e-173  0.000000e+00   0.000000e+00  0.000000e+00
29 29 -3.289785e-185  0.000000e+00   0.000000e+00  0.000000e+00
30 30 -4.906714e-198  0.000000e+00   0.000000e+00  0.000000e+00
31 31 -2.695250e-211  0.000000e+00   0.000000e+00  0.000000e+00
32 32 -5.452081e-225  0.000000e+00   0.000000e+00  0.000000e+00
33 33 -4.061186e-239  0.000000e+00   0.000000e+00  0.000000e+00
34 34 -1.113899e-253  0.000000e+00   0.000000e+00  0.000000e+00
35 35 -1.124911e-268  0.000000e+00   0.000000e+00  0.000000e+00
36 36 -4.182624e-284  0.000000e+00   0.000000e+00  0.000000e+00
37 37 -5.725571e-300  0.000000e+00   0.000000e+00  0.000000e+00
38 38 -2.885428e-316  0.000000e+00   0.000000e+00  0.000000e+00
39 39   0.000000e+00  0.000000e+00   0.000000e+00  0.000000e+00
40 40   0.000000e+00  0.000000e+00   0.000000e+00  0.000000e+00

`Phi(x)` Results:

get_phi = function(x) {
    data.frame(
        x = x,
        Phi_new = phi_new(x),
        Phi_old = phi_old(x),
        Phi_appr_new = phi_approx_new(x),
        Phi_appr_old = phi_approx_old(x)
    )
}

purrr::map_dfr(-40:9,get_phi)

     x       Phi_new       Phi_old  Phi_appr_new  Phi_appr_old
1  -40  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
2  -39  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
3  -38 2.885428e-316  0.000000e+00  0.000000e+00  0.000000e+00
4  -37 5.725571e-300 5.725571e-300  0.000000e+00  0.000000e+00
5  -36 4.182624e-284 4.182624e-284  0.000000e+00  0.000000e+00
6  -35 1.124911e-268 1.124911e-268  0.000000e+00  0.000000e+00
7  -34 1.113899e-253 1.113899e-253  0.000000e+00  0.000000e+00
8  -33 4.061186e-239 4.061186e-239  0.000000e+00  0.000000e+00
9  -32 5.452081e-225 5.452081e-225  0.000000e+00  0.000000e+00
10 -31 2.695250e-211 2.695250e-211  0.000000e+00  0.000000e+00
11 -30 4.906714e-198 4.906714e-198  0.000000e+00  0.000000e+00
12 -29 3.289785e-185 3.289785e-185  0.000000e+00  0.000000e+00
13 -28 8.123869e-173 8.123869e-173  0.000000e+00  0.000000e+00
14 -27 7.389481e-161 7.389481e-161  0.000000e+00  0.000000e+00
15 -26 2.476063e-149 2.476063e-149  0.000000e+00  0.000000e+00
16 -25 3.056697e-138 3.056697e-138  0.000000e+00  0.000000e+00
17 -24 1.390392e-127 1.390392e-127  0.000000e+00  0.000000e+00
18 -23 2.330637e-117 2.330637e-117  0.000000e+00  0.000000e+00
19 -22 1.439892e-107 1.439892e-107  0.000000e+00  0.000000e+00
20 -21  3.279278e-98  3.279278e-98 4.337000e-299 4.337000e-299
21 -20  2.753624e-89  2.753624e-89 9.394498e-260 9.394498e-260
22 -19  8.527224e-81  8.527224e-81 4.278579e-224 4.278579e-224
23 -18  9.740949e-73  9.740949e-73 6.256481e-192 6.256481e-192
24 -17  4.105996e-65  4.105996e-65 4.485680e-163 4.485680e-163
25 -16  6.388754e-58  6.388754e-58 2.408003e-137 2.408003e-137
26 -15  3.670966e-51  3.670966e-51 1.478016e-114 1.478016e-114
27 -14  7.793537e-45  7.793537e-45  1.584009e-94  1.584009e-94
28 -13  6.117164e-39  6.117164e-39  4.526424e-77  4.526424e-77
29 -12  1.776482e-33  1.776482e-33  5.266657e-62  5.266657e-62
30 -11  1.910660e-28  1.910660e-28  3.810306e-49  3.810306e-49
31 -10  7.619853e-24  7.619853e-24  2.617536e-38  2.617536e-38
32  -9  1.128588e-19  1.128588e-19  2.607333e-29  2.607333e-29
33  -8  6.220961e-16  6.220961e-16  5.750875e-22  5.750875e-22
34  -7  1.279813e-12  1.279813e-12  4.289120e-16  4.289120e-16
35  -6  9.865876e-10  9.865876e-10  1.651817e-11  1.651817e-11
36  -5  2.866516e-07  2.866516e-07  5.016240e-08  5.016240e-08
37  -4  3.167124e-05  3.167124e-05  1.834308e-05  1.834308e-05
38  -3  1.349898e-03  1.349898e-03  1.231955e-03  1.231955e-03
39  -2  2.275013e-02  2.275013e-02  2.276106e-02  2.276106e-02
40  -1  1.586553e-01  1.586553e-01  1.586697e-01  1.586697e-01
41   0  5.000000e-01  5.000000e-01  5.000000e-01  5.000000e-01
42   1  8.413447e-01  8.413447e-01  8.413303e-01  8.413303e-01
43   2  9.772499e-01  9.772499e-01  9.772389e-01  9.772389e-01
44   3  9.986501e-01  9.986501e-01  9.987680e-01  9.987680e-01
45   4  9.999683e-01  9.999683e-01  9.999817e-01  9.999817e-01
46   5  9.999997e-01  9.999997e-01  9.999999e-01  9.999999e-01
47   6  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00
48   7  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00
49   8  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00
50   9  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00

Also normal_lccdf should be updated too -- there might be a separate issue for that. Tnx for looking at this.

I'm not sure if there's a way to get log(Phi(x)) to replicate the precision of normal_lcdf(x|0,1). The precision is achieved through analytic approximations of the log function. Using log() severely reduces the precision and leads to faster overflow to 0. A log_Phi() function definitely makes sense though, as does updating normal_lccdf().

Also, I just happened to check the function reference doc and it still quotes the old values for underflow / overflow (-37.5 and 8.25 respectively) so this should probably be updated too.

Here is the docs issue :) stan-dev/docs#352

Navigation Menu

Repository navigation

Make `normal_lcdf(x | 0, 1)` equivalent to `Phi()` #2470

Testing

Model functions

`log(Phi(x))` Results:

`Phi(x)` Results:

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Make normal_lcdf(x | 0, 1) equivalent to Phi() #2470

Description

Activity

andrjohns commented on Apr 16, 2021

Testing

Model functions

log(Phi(x)) Results:

Phi(x) Results:

bbbales2 commented on Apr 16, 2021

PhilClemson commented on Apr 19, 2021

spinkney commented on Apr 19, 2021

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions

Make `normal_lcdf(x | 0, 1)` equivalent to `Phi()` #2470

`log(Phi(x))` Results:

`Phi(x)` Results: