Блог

Я обнаружил ошибку в реализации функции beta.ppf в scipy.stats. Это уже подтверждено и помечено как дефект в их системе исправления ошибок.

Однако в настоящее время мне нужно рассчитать доверительные интервалы бета-распределения, и поэтому мне нужна обратная бета-функция. Поскольку я не могу полагаться на текущую версию beta.ppf, мне нужна альтернатива для Python. Желательно, чтобы я не хотел внедрять функцию самостоятельно.

Кто-нибудь знает функцию, которая может заменить функцию beta.ppf из scipy.stats?

Итак, я создал обходной путь, используя библиотеку c boost, которая реализует версию бета-дистрибутива. Чтобы использовать его, вам необходимо установить библиотеку boost и скомпилировать ее. После этого вы можете использовать его, используя следующий код:

betainv.cpp

 #include <boost/python.hpp>
#include <boost/math/distributions/beta.hpp>

using namespace boost::python;

class betainvClass {
    public: double betainv(double p, double a, double b);
};

double betainvClass::betainv(double p, double a, double b) { 
    return boost::math::ibeta_inv(a, b, p);
}

// Expose classes and methods to Python
BOOST_PYTHON_MODULE(betainv) {
    class_<betainvClass> ("create_betainv_instance")
        .def("betainv", amp;betainvClass::betainv)
    ;
}
 

Makefile:

 TARGET = betainv
PYTHON_INC = /usr/local/Cellar/python/2.7.5/Frameworks/Python.framework/Versions/2.7/include/python2.7
PYTHON_LIB = /usr/local/Cellar/python/2.7.5/Frameworks/Python.framework/Versions/2.7/lib/python2.7/config
BOOST_INC = /usr/local/include
BOOST_LIB = /usr/local/lib

$(TARGET).so: $(TARGET).o
    g   -shared -Wl 
    $(TARGET).o -L$(BOOST_LIB) -lboost_python 
    -L$(PYTHON_LIB) -lpython2.7 
    -o $(TARGET).so

$(TARGET).o: $(TARGET).cpp
    g   -I$(PYTHON_INC) -I$(BOOST_INC) -c $(TARGET).cpp

clean:
    rm -f *.o *.a *.so *~ core
 

Файл примера Python:

 import betainv

beta = betainv.create_betainv_instance()
print "0.25, 0.0342, 170 -> "   str(beta.betainv(0.25, 0.0342, 170))
print "0.25, 0.0342, 171 -> "   str(beta.betainv(0.25, 0.0342, 171))
print "0.25, 0.0342, 172 -> "   str(beta.betainv(0.25, 0.0342, 172))
 

Кстати, еще один комментарий. Я провел тест скорости и запустил метод 1000 раз. Сначала с помощью scipy, а затем моей реализации boost. Это результаты в миллисекундах:

 required time scipy = 295.145019531
required time c     = 7.68383789062
 

ускорение реализации c происходит примерно в 42 раза быстрее.