用于极值的logit和inverse logit函数

Question

我需要logit和inverse logit函数logit(inv_logit(n)) == n.我使用numpy,这就是我所拥有的:

import numpy as np
def logit(p):
    return np.log(p) - np.log(1 - p)

def inv_logit(p):
    return np.exp(p) / (1 + np.exp(p))

以下是价值观:

print logit(inv_logit(2)) 
2.0 

print logit(inv_logit(10))
10.0 

print logit(inv_logit(20))
20.000000018 #well, pretty close

print logit(inv_logit(50))
Warning: divide by zero encountered in log
inf 

现在让我们测试负数

print logit(inv_logit(-10))
-10.0 
print logit(inv_logit(-20))
-20.0 
print logit(inv_logit(-200))
-200.0 
print logit(inv_logit(-500))
-500.0 
print logit(inv_logit(-2000))
Warning: divide by zero encountered in log
-inf 

所以我的问题是:实现这些功能的正确方法是什么,以便要求logit(inv_logit(n)) == n适用n于尽可能宽的范围(至少[-1e4; 1e4]？

而且(我确信这与第一个相关),为什么我的函数与正值相比更稳定？

Answer 1

要么使用

1. bigfloat包支持任意精度浮点运算.

2. SymPy 符号数学包.我将举两个例子:

首先,bigfloat:

http://packages.python.org/bigfloat/

这是一个简单的例子:

from bigfloat import *
def logit(p):
    with precision(100000):
        return log(p)- log(1 -BigFloat(p))

def inv_logit(p):
    with precision(100000):
        return exp(p) / (1 + exp(p))

int(round(logit(inv_logit(12422.0))))
# gives 12422
int(round(logit(inv_logit(-12422.0))))
# gives -12422

这真的很慢.您可能需要考虑重构问题并分析一些部分.像这样的案例在实际问题中很少见 - 我很好奇你正在处理什么样的问题.

安装示例:

wget http://pypi.python.org/packages/source/b/bigfloat/bigfloat-0.3.0a2.tar.gz
tar xvzf bigfloat-0.3.0a2.tar.gz 
cd bigfloat-0.3.0a2
as root:
python setup.py install

关于你的功能在负值下穿得更好的原因.考虑:

>>> float(inv_logit(-15))
3.059022269256247e-07

>>> float(inv_logit(15))
0.9999996940977731

在第一种情况下,浮点数很容易表示该值.移动小数点,以便不需要存储前导零:0.0000 .... 在第二种情况下,需要存储所有前导0.999,因此当稍后在logit()中执行1-p时,您需要所有额外的精度来获得精确的结果.

这是符号化的数学方式(明显更快!):

from sympy import *
def inv_logit(p):
    return exp(p) / (1 + exp(p))
def logit(p):
    return log(p)- log(1 -p)

x=Symbol('x')
expr=logit(inv_logit(x))
# expr is now:
# -log(1 - exp(x)/(1 + exp(x))) + log(exp(x)/(1 + exp(x)))
# rewrite it: (there are many other ways to do this. read the doc)
# you may want to make an expansion (of some suitable kind) instead.
expr=cancel(powsimp(expr)).expand()
# it is now 'x'

# just evaluate any expression like this:    
result=expr.subs(x,123.231)

# result is now an equation containing: 123.231
# to get the float: 
result.evalf()

Sympy在这里找到http://docs.sympy.org/.在ubuntu中,它是通过突触发现的.

Answer 2

有一种方法可以实现这些功能,使它们在各种值中都是稳定的,但它涉及根据参数区分的情况.

以inv_logit函数为例.你的公式"np.exp(p)/(1 + np.exp(p))"是正确的,但是对于大p会溢出.如果将分子和分母除以np.exp(p),则获得等效表达式

1. / (1. + np.exp(-p))

不同的是,这个不会溢出大正p.然而,对于p的大负值,它会溢出.因此,稳定的实施可以如下:

def inv_logit(p):
    if p > 0:
        return 1. / (1. + np.exp(-p))
    elif p <= 0:
        np.exp(p) / (1 + np.exp(p))
    else:
        raise ValueError

这是库LIBLINEAR(以及可能的其他库)中使用的策略.

Answer 3

现在，scipy有logit和expit（逆logit）函数，例如

>>> from scipy.special import logit, expit
>>> import numpy as np
>>> logit([0, 0.25, 0.5, 0.75, 1])
array([       -inf, -1.09861229,  0.        ,  1.09861229,         inf])
>>> expit([-np.inf, -1.5, 0, 1.5, np.inf])
array([ 0.        ,  0.18242552,  0.5       ,  0.81757448,  1.        ])