Блог

Я запускаю сэмплер MCMC, который требует вычисления гипергеометрической функции на каждом шаге с использованием scipy.special.hyp2f1() .

В определенных точках моей сетки (которые меня не волнуют) решения гипергеометрической функции довольно нестабильны, и SciPy выводит предупреждение:

 Warning! You should check the accuracy
  

Это довольно раздражает, и более 1000 выборок вполне могут замедлить мою работу.

Я безуспешно пытался использовать special.errprint(0) , а также отключил все предупреждения в Python, используя как warnings модуль, так и -W ignore флаг.

Нарушающая функция (вызываемая из другого файла) приведена ниже

 from numpy import pi, hypot, real, imag
import scipy.special as special

def deflection_angle(p, (x1, x2)):

    # Find the normalisation constant
    norm = (p.f * p.m * (p.r0 ** (t - 2.0)) / pi) ** (1.0 / t)

    # Define the complex plane
    z = x1   1j * x2

    # Define the radial coordinates
    r = hypot(x1, x2)

    # Truncate the radial coordinates
    r_ = r * (r < p.r0).astype('float')   p.r0 * (r >= p.r0).astype('float')

    # Calculate the radial part
    radial = (norm ** 2 / (p.f * z)) * ((norm / r_) ** (t - 2))

    # Calculate the angular part
    h1, h2, h3 = 0.5, 1.0 - t / 2.0, 2.0 - t / 2.0
    h4 = ((1 - p.f ** 2) / p.f ** 2) * (r_ / z) ** 2
    special.errprint(0)
    angular = special.hyp2f1(h1, h2, h3, h4)

    # Assemble the deflection angle
    alpha = (- radial * angular).conjugate()

    # Separate real and imaginary parts
    return real(alpha), imag(alpha)`
  

К сожалению, hyp2f1, как известно, сложно вычислить в некоторых нетривиальных областях пространства параметров. Многие реализации будут давать неточные или совершенно неправильные результаты. Scipy.special изо всех сил старается хотя бы отслеживать конвергенцию. Альтернативой может быть использование реализаций произвольной точности, например mpmath. Но это, безусловно, будет немного медленнее, поэтому пользователи MCMC остерегаются.

РЕДАКТИРОВАТЬ: Хорошо, похоже, это зависит от версии scipy. Я попробовал пример @wrwrwr в scipy 0.13.3, и он воспроизводит то, что вы видите: «Внимание! Вы должны проверить точность» печатается независимо от errprint статуса. Однако, делая то же самое с версией разработчика, я получаю

 In [12]: errprint(True)
Out[12]: 0

In [13]: hyp2f1(0.5, 2/3., 1.5, 0.09j 0.75j)
/home/br/virtualenvs/scipy_py27/bin/ipython:1: SpecialFunctionWarning: scipy.special/chyp2f1: loss of precision
  #!/home/br/virtualenvs/scipy_py27/bin/python
Out[13]: (0.93934867949609357 0.15593972567482395j)

In [14]: errprint(False)
Out[14]: 1

In [15]: hyp2f1(0.5, 2/3., 1.5, 0.09j 0.75j)
Out[15]: (0.93934867949609357 0.15593972567482395j)
  

Итак, по-видимому, это было исправлено в какой-то момент между 2013 и сейчас. Возможно, вы захотите обновить свою версию scipy.