Оцените производную отрицательного логарифмического правдоподобия 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.

Благодарю вас!


person njalex22    schedule 05.02.2021    source источник
comment
Прежде чем использовать lambdify, посмотрите на deriv. lambdify — это простой лексический транслятор, оставляющий такие вещи, как Derivative, Product и Dummy_39, которые не имеют смысла в numpy.   -  person hpaulj    schedule 05.02.2021


Ответы (2)


Во-первых, символический sympy и числовой numpy/scipy плохо сочетаются. Обычно рекомендуется разделять их. Кроме того, sympy не любит числа с плавающей запятой, поскольку они препятствуют точному символьному вычислению. Поэтому вместо np.pi предпочтительнее использовать символ sympy pi. Так же 0.5 нужно заменить на симпи дробь, например S.Half. (См. также подводные камни sympy. )

Производная от 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       
person JohanC    schedule 05.02.2021

Я знаю, что копирование и вставка выражений sympy не оптимальны, но мне нравится их видеть. Это дает читателям более четкое представление о том, что происходит. В противном случае им придется запускать код самостоятельно. Я могу сделать это, читая это на своем компьютере с подручным сеансом isympy, но не могу при использовании планшета или телефона.

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 раньше и даже не читал его документы, поэтому я не знаю, есть ли способ заставить его продолжить оценку. Если подумать, разве diff из f1(x)*f2(x)...*fn(x) не слишком длинное и грязное? Кое-что о правиле продукта для дифференциации.

Интересно, будет ли иметь значение установка 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 https://stackoverflow.com/a/66070692/901925:

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
person hpaulj    schedule 05.02.2021
comment
Спасибо. Кажется, это означает, что сумма не рассчитывается как сумма numpy. Таким образом, это может быть медленным для больших массивов. Я полагаю, нет никакого способа заставить lambdify создать массив numpy и суммировать? (Кстати, я только что узнал, что expand_log(..., force=True) создает сумму.) - person JohanC; 06.02.2021
comment
Я исследовал перемещение sum из lambdify, то есть получить функцию, которая может генерировать значения для каждого из test элементов и суммировать в конце. - person hpaulj; 06.02.2021