问题描述
背景:
我正在尝试实现一个执行
现在计算这个函数的逆:
x = sym.Symbol('x')逆 = sym.solve(sym.Eq(x, cdf), y)打印(反向)
输出:
[-LambertW((x - 1)*exp(-1)) - 1]
实际上,这只是给定 CDF 的负 y 的左分支:
sym.plot(inverse[0], (x, -0.5, 1))
问题:如何为给定 CDF 的正 y 获得正确的分支?
我的尝试:
指定
x
和y
为正:x = sym.Symbol('x', positive=True)y = sym.Symbol('y', positive=True)
这没有任何影响,即使对于第一个 CDF 图也是如此.
使 CDF 成为
Background:
I am trying to implement a function doing an inverse transform sampling. I use sympy for calculating CDF and getting its inverse function. While for some simple PDFs I get correct results, for a PDF which CDF's inverse function includes Lambert-W function, results are wrong.Example:
Consider following example CDF:import sympy as sym y = sym.Symbol('y') cdf = (-y - 1) * sym.exp(-y) + 1 # derived from `pdf = x * sym.exp(-x)` sym.plot(cdf, (y, -1, 5))
Now calculating inverse of this function:x = sym.Symbol('x') inverse = sym.solve(sym.Eq(x, cdf), y) print(inverse)
Output:
[-LambertW((x - 1)*exp(-1)) - 1]
This, in fact, is only a left branch of negative y's of a given CDF:
sym.plot(inverse[0], (x, -0.5, 1))
Question:How can I get the right branch for positive y's of a given CDF?
What I tried:
Specifying
x
andy
to be only positive:x = sym.Symbol('x', positive=True) y = sym.Symbol('y', positive=True)
This doesn't have any effect, even for the first CDF plot.
Making CDF a
Piecewise
function:cdf = sym.Piecewise((0, y < 0), ((-y - 1) * sym.exp(-y) + 1, True))
Again no effect. Strange thing here is that on another computer plotting this function gave a proper graph with zero for negative y's, but solving for a positive y's branch doesn't work anywhere. (Different versions? I also had to specify
adaptive=False
tosympy.plot
to make it work there.)Using
sympy.solveset
instead ofsympy.solve
:
This just gives a uselessConditionSet(y, Eq(x*exp(y) + y - exp(y) + 1, 0), Complexes(S.Reals x S.Reals, False))
as a result. Apparently,solveset
still doesn't know how to deal withLambertW
functions. From the docs:
Is it a bug or am I missing something? Is there any workaround to get the desired result?
解决方案The inverse produced by sympy is almost correct. The problem lies in the fact that the LambertW function has multiple branches over the domain (-1/e, 0). By default, it uses the upper branch, however for your problem you require the lower branch. The lower branch can be accessed by passing in a second argument to LambertW with a value of -1.
inverse = -sym.LambertW((x - 1)*sym.exp(-1), -1) - 1 sym.plot(inverse, (x, 0, 0.999))
Gives
这篇关于sympy.solve() 没有给出 LambertW 的解决方案之一的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!