Fondo:
Estoy tratando de implementar una función haciendo una transformación inversa de muestreo . Uso sympy para calcular CDF y obtener su función inversa. Mientras que para algunos PDF simples obtengo resultados correctos, para un PDF cuya función inversa de CDF incluye la función de Lambert-W , los resultados son incorrectos.
Ejemplo:
Considere el siguiente ejemplo de 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))
Ahora calculando el inverso de esta función:
x = sym.Symbol('x') inverse = sym.solve(sym.Eq(x, cdf), y) print(inverse)
Salida:
[-LambertW((x - 1)*exp(-1)) - 1]
Esto, de hecho, es solo una twig izquierda de las y negativas de un CDF dado:
sym.plot(inverse[0], (x, -0.5, 1))
Pregunta: ¿Cómo puedo obtener la twig correcta para las y positivas de un CDF dado?
Lo que intenté:
Especificando x
es solo positivo:
x = sym.Symbol('x', positive=True) y = sym.Symbol('y', positive=True)
Esto no tiene ningún efecto, incluso para la primera gráfica de la FCD.
Haciendo de CDF una función a Piecewise
:
cdf = sym.Piecewise((0, y < 0), ((-y - 1) * sym.exp(-y) + 1, True))
De nuevo sin efecto. Lo extraño aquí es que en otra computadora, al trazar esta función se obtuvo un gráfico adecuado con cero para las y negativas, pero la resolución de una twig positiva de y no funciona en ninguna parte. (¿Versiones diferentes? También tuve que especificar adaptive=False
to sympy.plot
para que funcione allí).
Usando sympy.solveset
lugar de sympy.solve
:
Esto solo da un ConditionSet(y, Eq(x*exp(y) + y - exp(y) + 1, 0), Complexes(S.Reals x S.Reals, False))
inútil ConditionSet(y, Eq(x*exp(y) + y - exp(y) + 1, 0), Complexes(S.Reals x S.Reals, False))
como resultado. Aparentemente, solveset
aún no sabe cómo lidiar con LambertW
funciones de LambertW
. De los documentos:
Cuando los casos no se resuelven o solo se pueden resolver de manera incompleta, se utiliza un
ConditionSet
y actúa como un objeto de Solveset sin evaluar. Todavía hay algunas cosas quesolveset
no puede hacer, lo que puedesolve
el antiguo, como resolver ecuaciones de tipo multivariable y LambertW no lineales.
¿Es un error o me falta algo? ¿Hay alguna solución para obtener el resultado deseado?
La inversa producida por sympy es casi correcta. El problema radica en el hecho de que la función LambertW tiene varias twigs sobre el dominio (-1 / e, 0). De forma predeterminada, utiliza la twig superior, sin embargo, para su problema, necesita la twig inferior. Se puede acceder a la twig inferior pasando un segundo argumento a LambertW con un valor de -1.
inverse = -sym.LambertW((x - 1)*sym.exp(-1), -1) - 1 sym.plot(inverse, (x, 0, 0.999))
Da