评估负对数似然的导数 Python
Evaluate Derivative of Negative Log Likelihood Python
我正在尝试计算 python 中负对数似然函数的导数。我正在使用 sympy 来计算导数，但是，当我尝试对其求值时收到错误消息。下面的示例代码尝试为对数正态函数计算此值。
import numpy as np
import sympy as sym
from sympy import Product, Function, oo, IndexedBase, diff, Eq, symbols, log, exp
from scipy.stats import lognorm
np.random.seed(seed=111)
test = lognorm.rvs(s=1,loc=2,scale=1,size=1000)
x = IndexedBase('x')
i = symbols('i', positive=True)
n = symbols('n', positive=True)
mx = symbols('mx', positive=True)
sx = symbols('sx', positive=True)
pdf = 1 / (x[i] * sx * sqrt(2 * np.pi)) * exp(-0.5 * ((log(x[i]-mx))**2/(sx)**2))
Log_LL = -log(Product(pdf, (i, 1, n)))
deriv = diff(Log_LL, mx)
fx = lambdify([x,mx,sx,n],deriv)
fx(test,2,1,len(test))
当我计算这个公式时，我收到以下错误：
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-139-452f00dee1f1> in <module>
1 fx = lambdify([x,mx,sx,n],deriv)
----> 2 fx(test,2,1,len(test))
<lambdifygenerated-14> in _lambdifygenerated(Dummy_39, mx, sx, n)
3 # Derivative
4 # Product
----> 5 -Derivative(Product(0.398942280401433*exp(-0.5*log(-mx + Dummy_39[i])**2/sx**2)/(sx*Dummy_39[i]), (i, 1, n)), mx)/Product(0.398942280401433*exp(-0.5*log(-mx + Dummy_39[i])**2/sx**2)/(sx*Dummy_39[i]), (i, 1, n)))
NameError: name 'Derivative' is not defined
我知道表达式“deriv”包含导数运算符，但我只是计算单个变量的导数，所以我相信 sympy 应该能够处理这个问题。
我是 运行 sympy 1.7.1、numpy 1.19.2 和 scipy 1.5.2
谢谢！
我知道 sympy 表达式的复制粘贴显示不是最佳的，但我喜欢看到它们。它让读者更清楚地了解正在发生的事情。否则他们必须自己 运行 代码。我可以在电脑上阅读这篇文章时使用 isympy 方便的会话，但在使用平板电脑或 phone.
时不能
In [164]: pdf
Out[164]:
2
-0.5⋅log (-mx + x[i])
──────────────────────
2
sx
0.398942280401433⋅ℯ
─────────────────────────────────────────
sx⋅x[i]
In [165]: Log_LL
Out[165]:
⎛ n ⎞
⎜─┬────────────┬─ ⎟
⎜ │ │ 2 ⎟
⎜ │ │ -0.5⋅log (-mx + x[i]) ⎟
⎜ │ │ ──────────────────────⎟
⎜ │ │ 2 ⎟
-log⎜ │ │ sx ⎟
⎜ │ │ 0.398942280401433⋅ℯ ⎟
⎜ │ │ ─────────────────────────────────────────⎟
⎜ │ │ sx⋅x[i] ⎟
⎜ │ │ ⎟
⎝ i = 1 ⎠
所以 log 包裹了 Product（用大 pi 符号显示），但不尝试对内部表达式做任何事情。
In [166]: deriv = diff(Log_LL, mx)
In [167]: deriv
Out[167]:
⎛ n ⎞
⎜─┬────────────┬─ ⎟
⎜ │ │ 2 ⎟
⎜ │ │ -0.5⋅log (-mx + x[i]) ⎟
⎜ │ │ ──────────────────────⎟
∂ ⎜ │ │ 2 ⎟
-───⎜ │ │ sx ⎟
∂mx⎜ │ │ 0.398942280401433⋅ℯ ⎟
⎜ │ │ ─────────────────────────────────────────⎟
⎜ │ │ sx⋅x[i] ⎟
⎜ │ │ ⎟
⎝ i = 1 ⎠
─────────────────────────────────────────────────────────────────
n
─┬────────────┬─
│ │ 2
│ │ -0.5⋅log (-mx + x[i])
│ │ ──────────────────────
│ │ 2
│ │ sx
│ │ 0.398942280401433⋅ℯ
│ │ ─────────────────────────────────────────
│ │ sx⋅x[i]
│ │
i = 1
diff 执行 df(x)/dx/x 对数微分，但不会尝试作用于 Product。
我以前没有用过Product，甚至没有读过它的文档，所以我不知道是否有办法强制进一步评估。想想看，f1(x)*f2(x)...*fn(x)的diff是不是又长又乱？关于差异化的产品规则。
我想知道将 n 设置为一个小数字（例如 3）而不是一个变量是否会有所不同。
无论如何，通过 lambdify 传递它是没有意义的。正如错误消息中的代码所示，lambdify 不会评估或清理 sympy，它只是进行简单的词法替换。
-Derivative(Product(0.398942280401433*exp(-0.5*log(-mx + Dummy_39[i])**2/sx**2) /
(sx*Dummy_39[i]), (i, 1, n)), mx)/
Product(0.398942280401433*exp(-0.5*log(-mx + Dummy_39[i])**2/sx**2) /
(sx*Dummy_39[i]), (i, 1, n)))
===
经过 JohanC 的更改 :
In [182]: deriv
Out[182]:
n - 1
____
╲
╲ log(-mx + x[i])
╲ ────────────────
- ╱ 2
╱ sx ⋅(-mx + x[i])
╱
‾‾‾‾
i = 0
In [183]: fx = lambdify([x, mx, sx, n], deriv)
...:
In [184]: fx??
Signature: fx(Dummy_2642, mx, sx, n)
Docstring:
Created with lambdify. Signature:
func(x, mx, sx, n)
Expression:
-Sum(log(-mx + x[i])/(sx**2*(-mx + x[i])), (i, 0, n - 1))
Source code:
def _lambdifygenerated(Dummy_2642, mx, sx, n):
return (-(builtins.sum(log(-mx + Dummy_2642[i])/(sx**2*(-mx + Dummy_2642[i])) for i in range(0, n - 1+1))))
===
根据 JohanC 的评论，该函数可以（应该？）重写为
def foo(x, mx, sx, n):
return -np.sum(np.log(-mx + x)/(sx**2*(-mx + x)))
对于小样本：
In [206]: test
Out[206]:
array([2.32179572, 3.46861414, 6.46627075, 2.70090544, 2.45496557,
2.63163795, 2.94254768, 2.7017532 , 2.47925472, 2.30607048])
In [207]: fx(test,2,1,len(test))
Out[207]: 11.862478799879577
In [208]: foo(test,2,1,len(test))
Out[208]: 11.862478799879577
显然 lambdify 没有足够的 numpy“知识”来用“矢量化”解决方案替换 sum('comprehension')。
但是我们可以利用我们自己的知识为个体x[i]生成值，这里用一个通用的y符号代替：
In [224]: pdfy = 1 / (y * sx * sqrt(2 * pi)) * exp(-S.Half * ((log(y - mx)) ** 2 / (sx) ** 2))
In [225]: f1 = diff(log(pdfy), mx)
In [226]: fx1 = lambdify([y, mx,sx], f1)
In [227]: fx1?
Signature: fx1(y, mx, sx)
Docstring:
Created with lambdify. Signature:
func(y, mx, sx)
Expression:
log(-mx + y)/(sx**2*(-mx + y))
Source code:
def _lambdifygenerated(y, mx, sx):
return (log(-mx + y)/(sx**2*(-mx + y)))
In [228]: fx1(test, 2, 1)
Out[228]:
array([-3.52347235, 0.26168834, 0.33507905, -0.50703316, -1.73097394,
-0.72737698, -0.06277536, -0.50469809, -1.53472262, -3.86819368])
In [229]: -fx1(test, 2, 1).sum()
Out[229]: 11.862478799879577
首先，符号 sympy 和数字 numpy/scipy 不能很好地混合。通常建议将它们分开。此外，sympy 不喜欢浮点数，因为它们会阻止精确的符号计算。因此，首选 sympy 的符号 pi 而不是 np.pi。此外，0.5 应替换为 sympy 分数，例如 S.Half。 （另请参阅 sympy's gotchas。）
Sympy 的导数似乎无法应对Product。我们可以尝试用日志的总和来代替乘积的日志。 expand_log(..., force=True) 可以帮助进行该转换（force=True 当 sympy 不确定表达式是否肯定为正时，大概 x[i] 可能很复杂）。
转为numpy时，numpy不喜欢从1开始索引，可以通过从0开始索引解决
from sympy import Product, Sum, IndexedBase, diff, symbols, log, exp, lambdify, sqrt, pi, S, expand_log
x = IndexedBase('x')
i = symbols('i', positive=True)
n = symbols('n', positive=True)
mx = symbols('mx', positive=True)
sx = symbols('sx', positive=True)
pdf = 1 / (x[i] * sx * sqrt(2 * pi)) * exp(-S.Half * ((log(x[i] - mx)) ** 2 / (sx) ** 2))
# Log_LL = -Sum(log(pdf), (i, 0, n - 1))
Log_LL = -log(Product(pdf, (i, 0, n-1)))
Log_LL = expand_log(Log_LL, force=True)
deriv = diff(Log_LL, mx)
fx = lambdify([x, mx, sx, n], deriv)
from scipy.stats import lognorm
import numpy as np
np.random.seed(seed=111)
test = lognorm.rvs(s=1, loc=2, scale=1, size=1000)
fx(test, 2, 1, len(test))
这样，sympy 设法以符号方式计算导数：
n - 1
____
╲
╲ log(-mx + x[i])
╲ ────────────────
- ╱ 2
╱ sx ⋅(-mx + x[i])
╱
‾‾‾‾
i = 0
我正在尝试计算 python 中负对数似然函数的导数。我正在使用 sympy 来计算导数，但是，当我尝试对其求值时收到错误消息。下面的示例代码尝试为对数正态函数计算此值。
import numpy as np
import sympy as sym
from sympy import Product, Function, oo, IndexedBase, diff, Eq, symbols, log, exp
from scipy.stats import lognorm
np.random.seed(seed=111)
test = lognorm.rvs(s=1,loc=2,scale=1,size=1000)
x = IndexedBase('x')
i = symbols('i', positive=True)
n = symbols('n', positive=True)
mx = symbols('mx', positive=True)
sx = symbols('sx', positive=True)
pdf = 1 / (x[i] * sx * sqrt(2 * np.pi)) * exp(-0.5 * ((log(x[i]-mx))**2/(sx)**2))
Log_LL = -log(Product(pdf, (i, 1, n)))
deriv = diff(Log_LL, mx)
fx = lambdify([x,mx,sx,n],deriv)
fx(test,2,1,len(test))
当我计算这个公式时，我收到以下错误：
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-139-452f00dee1f1> in <module>
1 fx = lambdify([x,mx,sx,n],deriv)
----> 2 fx(test,2,1,len(test))
<lambdifygenerated-14> in _lambdifygenerated(Dummy_39, mx, sx, n)
3 # Derivative
4 # Product
----> 5 -Derivative(Product(0.398942280401433*exp(-0.5*log(-mx + Dummy_39[i])**2/sx**2)/(sx*Dummy_39[i]), (i, 1, n)), mx)/Product(0.398942280401433*exp(-0.5*log(-mx + Dummy_39[i])**2/sx**2)/(sx*Dummy_39[i]), (i, 1, n)))
NameError: name 'Derivative' is not defined
我知道表达式“deriv”包含导数运算符，但我只是计算单个变量的导数，所以我相信 sympy 应该能够处理这个问题。
我是 运行 sympy 1.7.1、numpy 1.19.2 和 scipy 1.5.2
谢谢！
我知道 sympy 表达式的复制粘贴显示不是最佳的，但我喜欢看到它们。它让读者更清楚地了解正在发生的事情。否则他们必须自己 运行 代码。我可以在电脑上阅读这篇文章时使用 isympy 方便的会话，但在使用平板电脑或 phone.
In [164]: pdf
Out[164]:
2
-0.5⋅log (-mx + x[i])
──────────────────────
2
sx
0.398942280401433⋅ℯ
─────────────────────────────────────────
sx⋅x[i]
In [165]: Log_LL
Out[165]:
⎛ n ⎞
⎜─┬────────────┬─ ⎟
⎜ │ │ 2 ⎟
⎜ │ │ -0.5⋅log (-mx + x[i]) ⎟
⎜ │ │ ──────────────────────⎟
⎜ │ │ 2 ⎟
-log⎜ │ │ sx ⎟
⎜ │ │ 0.398942280401433⋅ℯ ⎟
⎜ │ │ ─────────────────────────────────────────⎟
⎜ │ │ sx⋅x[i] ⎟
⎜ │ │ ⎟
⎝ i = 1 ⎠
所以 log 包裹了 Product（用大 pi 符号显示），但不尝试对内部表达式做任何事情。
In [166]: deriv = diff(Log_LL, mx)
In [167]: deriv
Out[167]:
⎛ n ⎞
⎜─┬────────────┬─ ⎟
⎜ │ │ 2 ⎟
⎜ │ │ -0.5⋅log (-mx + x[i]) ⎟
⎜ │ │ ──────────────────────⎟
∂ ⎜ │ │ 2 ⎟
-───⎜ │ │ sx ⎟
∂mx⎜ │ │ 0.398942280401433⋅ℯ ⎟
⎜ │ │ ─────────────────────────────────────────⎟
⎜ │ │ sx⋅x[i] ⎟
⎜ │ │ ⎟
⎝ i = 1 ⎠
─────────────────────────────────────────────────────────────────
n
─┬────────────┬─
│ │ 2
│ │ -0.5⋅log (-mx + x[i])
│ │ ──────────────────────
│ │ 2
│ │ sx
│ │ 0.398942280401433⋅ℯ
│ │ ─────────────────────────────────────────
│ │ sx⋅x[i]
│ │
i = 1
diff 执行 df(x)/dx/x 对数微分，但不会尝试作用于 Product。
我以前没有用过Product，甚至没有读过它的文档，所以我不知道是否有办法强制进一步评估。想想看，f1(x)*f2(x)...*fn(x)的diff是不是又长又乱？关于差异化的产品规则。
我想知道将 n 设置为一个小数字（例如 3）而不是一个变量是否会有所不同。
无论如何，通过 lambdify 传递它是没有意义的。正如错误消息中的代码所示，lambdify 不会评估或清理 sympy，它只是进行简单的词法替换。
-Derivative(Product(0.398942280401433*exp(-0.5*log(-mx + Dummy_39[i])**2/sx**2) /
(sx*Dummy_39[i]), (i, 1, n)), mx)/
Product(0.398942280401433*exp(-0.5*log(-mx + Dummy_39[i])**2/sx**2) /
(sx*Dummy_39[i]), (i, 1, n)))
===
经过 JohanC 的更改
In [182]: deriv
Out[182]:
n - 1
____
╲
╲ log(-mx + x[i])
╲ ────────────────
- ╱ 2
╱ sx ⋅(-mx + x[i])
╱
‾‾‾‾
i = 0
In [183]: fx = lambdify([x, mx, sx, n], deriv)
...:
In [184]: fx??
Signature: fx(Dummy_2642, mx, sx, n)
Docstring:
Created with lambdify. Signature:
func(x, mx, sx, n)
Expression:
-Sum(log(-mx + x[i])/(sx**2*(-mx + x[i])), (i, 0, n - 1))
Source code:
def _lambdifygenerated(Dummy_2642, mx, sx, n):
return (-(builtins.sum(log(-mx + Dummy_2642[i])/(sx**2*(-mx + Dummy_2642[i])) for i in range(0, n - 1+1))))
===
根据 JohanC 的评论，该函数可以（应该？）重写为
def foo(x, mx, sx, n):
return -np.sum(np.log(-mx + x)/(sx**2*(-mx + x)))
对于小样本：
In [206]: test
Out[206]:
array([2.32179572, 3.46861414, 6.46627075, 2.70090544, 2.45496557,
2.63163795, 2.94254768, 2.7017532 , 2.47925472, 2.30607048])
In [207]: fx(test,2,1,len(test))
Out[207]: 11.862478799879577
In [208]: foo(test,2,1,len(test))
Out[208]: 11.862478799879577
显然 lambdify 没有足够的 numpy“知识”来用“矢量化”解决方案替换 sum('comprehension')。
但是我们可以利用我们自己的知识为个体x[i]生成值，这里用一个通用的y符号代替：
In [224]: pdfy = 1 / (y * sx * sqrt(2 * pi)) * exp(-S.Half * ((log(y - mx)) ** 2 / (sx) ** 2))
In [225]: f1 = diff(log(pdfy), mx)
In [226]: fx1 = lambdify([y, mx,sx], f1)
In [227]: fx1?
Signature: fx1(y, mx, sx)
Docstring:
Created with lambdify. Signature:
func(y, mx, sx)
Expression:
log(-mx + y)/(sx**2*(-mx + y))
Source code:
def _lambdifygenerated(y, mx, sx):
return (log(-mx + y)/(sx**2*(-mx + y)))
In [228]: fx1(test, 2, 1)
Out[228]:
array([-3.52347235, 0.26168834, 0.33507905, -0.50703316, -1.73097394,
-0.72737698, -0.06277536, -0.50469809, -1.53472262, -3.86819368])
In [229]: -fx1(test, 2, 1).sum()
Out[229]: 11.862478799879577
首先，符号 sympy 和数字 numpy/scipy 不能很好地混合。通常建议将它们分开。此外，sympy 不喜欢浮点数，因为它们会阻止精确的符号计算。因此，首选 sympy 的符号 pi 而不是 np.pi。此外，0.5 应替换为 sympy 分数，例如 S.Half。 （另请参阅 sympy's gotchas。）
Sympy 的导数似乎无法应对Product。我们可以尝试用日志的总和来代替乘积的日志。 expand_log(..., force=True) 可以帮助进行该转换（force=True 当 sympy 不确定表达式是否肯定为正时，大概 x[i] 可能很复杂）。
转为numpy时，numpy不喜欢从1开始索引，可以通过从0开始索引解决
from sympy import Product, Sum, IndexedBase, diff, symbols, log, exp, lambdify, sqrt, pi, S, expand_log
x = IndexedBase('x')
i = symbols('i', positive=True)
n = symbols('n', positive=True)
mx = symbols('mx', positive=True)
sx = symbols('sx', positive=True)
pdf = 1 / (x[i] * sx * sqrt(2 * pi)) * exp(-S.Half * ((log(x[i] - mx)) ** 2 / (sx) ** 2))
# Log_LL = -Sum(log(pdf), (i, 0, n - 1))
Log_LL = -log(Product(pdf, (i, 0, n-1)))
Log_LL = expand_log(Log_LL, force=True)
deriv = diff(Log_LL, mx)
fx = lambdify([x, mx, sx, n], deriv)
from scipy.stats import lognorm
import numpy as np
np.random.seed(seed=111)
test = lognorm.rvs(s=1, loc=2, scale=1, size=1000)
fx(test, 2, 1, len(test))
这样，sympy 设法以符号方式计算导数：
n - 1
____
╲
╲ log(-mx + x[i])
╲ ────────────────
- ╱ 2
╱ sx ⋅(-mx + x[i])
╱
‾‾‾‾
i = 0