Adding log(ive) Modified Bessel Function of the first kind

[sethtroisi]

As a follow up to #12244 I had fun adding a function to calculate log(ive).

@ev-br mentioned this is already done in _hyp0f1.pxd#L51
This points to DLMF 10.41.10, I extended _hyp0f1 to k=5 using A144617, which reduces max error.

I'd love some advice from @WarrenWeckesser about

• where this code belongs
• if it needs to support complex
• how to make sure we can call it from stats.ncx2

---

using timeit in colab, special.ive takes ~1.5 µs, new cython code takes ~500 ns, pure python takes 18 µs

Error is still large when v and z are small. It probably makes sense to call np.log(special.ive(v, z)) when v and z are small and only use this when v and z are large.

Click to toggle implementation

import scipy.special
import numpy as np

def test_ive_asym(v, z):
    r"""Asymptotic expansion for I_{v}(z)
    for real $z > 0$ and $v\to +\infty$.
    Based off DLMF 10.41
    """
    # DLMF 10.41.3 calculates I_v(v*z)
    x = z / v

    # DLMF 10.41.8
    p = np.sqrt(1.0 + x*x)

    # DLMF 10.41.7
    eta = p + np.log(x) - scipy.special.log1p(p)

    # e^(vn) / ((2piv)^(1/2)*p^(1/2))
    log_res = v*eta - (np.log(2.0*np.pi*v)/2 + np.log(p)/2)

    # large-v asymptotic correction, DLMF 10.41.10
    pp = 1.0/p
    p2 = pp*pp
    p4 = p2*p2
    p6 = p4*p2
    p8 = p4*p4
    p10 = p6*p4

    v2 = v*v
    v3 = v2*v
    v4 = v2*v2
    v5 = v3*v2

    # u1,u2,u3 from DLMF, higher terms from oeis.org/A144617 and oeis.org/A144618
    u1 = (3.0 - 5.0*p2) * pp / 24.0
    u2 = (81.0 - 462.0*p2 + 385.0*p4) * p2 / 1152.0
    u3 = (30375 - 369603*p2 + 765765*p4 - 425425*p6) * pp * p2 / 414720.0
    u4 = (4465125 - 94121676*p2 + 349922430*p4 - 446185740*p6 + 185910725*p8) * p4 / 39813120.0
    u5 = (1519035525.0 - 49286948607.0*p2 + 284499769554.0*p4   - 614135872350.0*p6 + 566098157625.0*p8 - 188699385875.0*p10) * pp * p4 / 6688604160.0
    u_corr_i = 1.0 + u1/v + u2/v2 + u3/v3 + u4/v4 + u5/v5

    # -z for ive
    result = log_res + np.log(u_corr_i) - abs(z)
    return result


for i, (v, z) in enumerate([
    (3, 1), (3, 5), (3, 10),
    (3, 50), (3, 100), (3, 200),
    (10, 1), (10, 5), (10, 10),
    (50, 1), (50, 5), (50, 10),
    (150, 1), (150, 5), (150, 10),
]):
    with np.errstate(all='ignore'):
        ive = scipy.special.ive(v, z)
        log_ive = np.log(ive)
        log_ive2 = test_ive_asym(v, z)
        err = (log_ive2 - log_ive) / log_ive
        print(f"  ive({v:3}, {z:2}) = {ive:9.2e}, log = {log_ive:.2f} vs {log_ive2:.2f} err: {err:.1e}")
        if i % 3 == 2: print()

Click to toggle results

  ive(  3,  1) =  8.16e-03, log = -4.81 vs -4.81 err: 4.5e-06
  ive(  3,  5) =  6.96e-02, log = -2.66 vs -2.66 err: 3.5e-06
  ive(  3, 10) =  7.98e-02, log = -2.53 vs -2.53 err: -2.0e-07

  ive(  3, 50) =  5.16e-02, log = -2.96 vs -2.96 err: 1.1e-11
  ive(  3, 100) =  3.82e-02, log = -3.27 vs -3.27 err: 1.7e-13
  ive(  3, 200) =  2.76e-02, log = -3.59 vs -3.59 err: 7.1e-15

  ive( 10,  1) =  1.01e-10, log = -23.01 vs -23.01 err: -4.9e-11
  ive( 10,  5) =  3.09e-05, log = -10.39 vs -10.39 err: -6.6e-10
  ive( 10, 10) =  9.94e-04, log = -6.91 vs -6.91 err: 8.2e-10

  ive( 50,  1) =  1.08e-80, log = -184.13 vs -184.13 err: -0.0e+00
  ive( 50,  5) =  1.98e-47, log = -107.54 vs -107.54 err: -5.3e-16
  ive( 50, 10) =  2.16e-34, log = -77.52 vs -77.52 err: 3.8e-15

  ive(150,  1) =  0.00e+00, log = -inf vs -709.99 err: nan
  ive(150,  5) = 6.03e-206, log = -472.54 vs -472.54 err: -0.0e+00
  ive(150, 10) = 6.57e-163, log = -373.44 vs -373.44 err: 1.5e-16

[sethtroisi]

In colab format

[jvavrek]

This could help with the underflow in the log term of #11777

[dschmitz89]

log(ive) would also be very useful for the von Mises distribution (for example #17624).

If anyone wants to work on this, it might be worth looking into tensorflow which has an implementation: https://github.com/tensorflow/probability/blob/0759c57eb0306a5bb0fcc3d105bbcab9943092a5/tensorflow_probability/python/math/bessel.py#L879 .