Проблема определения специальных функций в Python против Mathematica - PullRequest
Новая
       89

Проблема определения специальных функций в Python против Mathematica

2 голосов
/ 22 марта 2019

У меня есть код Mathematica, который вычисляет 95% доверительные интервалы накопительной функции распределения (CDF), полученной из определенной функции распределения вероятности (PDF).PDF-файл ужасен, так как содержит гипергеометрическую функцию 2F1, и мне нужно вычислить 2-сигма-погрешности набора данных из 15 значений.

Я хочу перевести этот код на Python, но я получаю очень существенное расхождение во второй половине значений.

Код Mathematica

results нижеи верхний уровень достоверности 2-сигма для значений в xdata.То есть xdata всегда должно находиться между двумя соответствующими results значениями. 

navs  = {10, 10, 18, 30, 52, 87, 147, 245, 410, 684, 1141, 1903, 3173,  5290, 8816};
freqs = {0.00002, 0.00004, 0.0000666667, 0.000111111, 0.000185185,   0.000308642, 0.000514403, 0.000857339, 0.00142893, 0.00238166,   0.00396944, 0.00661594, 0.0165426, 0.0220568, 0.027571}
xdata = {0.578064980346793,   0.030812200935204,   0.316777979844816,  
         0.353718150091612,   0.287659600326548,   0.269254388840293,  
         0.16545714457921,   0.138759871084825,    0.0602382519940077,   
         0.10120771961,  0.065311134782518,    0.105235790998594,   
         0.124642033979457,   0.0271909963701794,  0.0686653810421847};
data = MapThread[{#1, #2, #3} &, {navs, freqs, xdata}]

post[x_, n_, y_] = 
     (n - 1) (1 - x)^n (1 - y)^(n - 2) Hypergeometric2F1[n, n, 1, x*y]

integral = Map[(values = #; mesh = Subdivide[0, 1, 1000]; 
     Interpolation[
      DeleteDuplicates[{Map[
            SetPrecision[post[#, values[[1]], values[[3]]^2], 100] &, 
            mesh] // (Accumulate[#] - #/2 - #[[1]]/
               2) & // #/#[[-1]] &, 
         mesh}\[Transpose], (#1[[1]] == #2[[1]] &)], 
      InterpolationOrder -> 1]) &, data];

results = 
 MapThread[{Sqrt[#1[.025]], Sqrt[#1[0.975]]} &, {integral, data}]

{{0.207919, 0.776508}, {0.0481485, 0.535278}, {0.0834002, 0.574447}, 
{0.137742, 0.551035}, {0.121376, 0.455097}, {0.136889, 0.403306}, 
{0.0674029, 0.279408}, {0.0612534, 0.228762}, {0.0158357, 0.134521}, 
{0.0525374, 0.156055}, {0.0270589, 0.108861}, {0.0740978, 0.137691}, 
{0.100498, 0.149646}, {0.00741129, 0.0525161}, {0.0507748, 0.0850961}}

Код Python

Вот мой перевод: results - это то же количество, что и раньше, усеченноедо 7-й цифры для повышения читабельности.

Получаемые мной значения results начинают расходиться с 7-й пары значений, а последние четыре точки xdata do not попадают между двумя соответствующими results значениями.

import numpy as np
from scipy.integrate import cumtrapz
from scipy.interpolate import interp1d
from mpmath import *

mesh = list(np.linspace(0,1,1000));

navs = [10, 10, 18, 30, 52, 87, 147, 245, 410, 684, 1141, 1903, 3173, 5290, 8816]
freqs = [0.00002, 0.00004, 0.0000666667, 0.000111111, 0.000185185, 0.000308642, 0.000514403, 0.000857339, 0.00142893, 0.00238166, 0.00396944, 0.00661594, 0.0165426, 0.0220568, 0.027571]
xdata = [0.578064980346793, 0.030812200935204, 0.316777979844816, 
0.353718150091612,0.287659600326548, 0.269254388840293,
0.16545714457921, 0.138759871084825, 0.0602382519940077, 
0.10120771961, 0.065311134782518, 0.105235790998594, 
0.124642033979457, 0.0271909963701794, 0.0686653810421847]

def post(x,n,y):
    post = (n-1)*((1-x)**n)*((1-y)**(n-2))*hyp2f1(n,n,1,x*y)
    return post

# setting the numeric precision to 100 as in Mathematica
# trying to get the most precise hypergeometric function values
mp.dps = 100
mp.pretty = True

results = []

for i in range(len(navs)):
    postprob = [];
    for j in range(len(mesh)):    
        posterior = post(mesh[j], navs[i], xdata[i]**2)
        postprob.append(posterior)
# calculate the norm of the pdf for integration
    norm = np.trapz(np.array(postprob),mesh);
# integrate pdf/norm to obtain cdf
    integrate = list(np.unique(cumtrapz(np.array(postprob)/norm, mesh, initial=0)));
    mesh2 = list(np.linspace(0,1,len(integrate)));
# interpolate inverse cdf to obtain the 2sigma quantiles
    icdf = interp1d(integrate, mesh2, bounds_error=False, fill_value='extrapolate');
    results.append(list(np.sqrt(icdf([0.025, 0.975]))))

results

[[0.2079198, 0.7765088], [0.0481485, 0.5352773], [0.0834, 0.5744489],
 [0.1377413, 0.5510352], [0.1218029, 0.4566994], [0.1399324, 0.4122767],
 [0.0733743, 0.3041607], [0.0739691, 0.2762597], [0.0230135, 0.1954886],
 [0.0871462, 0.2588804], [0.05637, 0.2268962],   [0.1731199, 0.3217401],
 [0.2665897, 0.3969059], [0.0315915, 0.2238736], [0.2224567, 0.3728803]]

Благодаря комментариям к этому вопросу я обнаружил, что:

  • Гипергеометрическая функция дает разные результаты на двух языках.С теми же входными значениями я получаю это: В Mathematica Hypergeometric2F1 дает мне в результате 1.0588267, в то время как в Python mpmath.hyp2f1 дает 1.0588866.Это самая вторая точка сетки, и разница в пятом десятичном знаке.

Есть ли где-то лучшее определение этой специальной функции, которую я не смог найти?

  • Я до сих пор не знаю, происходит ли это только из-за гипергеометрической функции или также из-за метода интегрирования, но это определенно отправная точка.

(Я довольно новичок вPython, может быть, код немного наивен)

...