python - 极值的 logit 和反 logit 函数

Question

我需要 logit 和反 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

在第一种情况下， float 很容易表示这个值。移动小数点以便前导零: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 中，它是通过突触发现的。