Почему в SymPy суммирование 2-х случайных величин дает непостижимый результат?

Question

С учетом суммирования 2 гауссовых переменных в SymPy:

from sympy import *
from sympy.stats import *

init_printing()

a = Normal('a', 0, Symbol('P', real=true))
b = Normal('b', 0, Symbol('Q', real=true))

ss = a + b

pprint(simplify(density(ss)))

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

    ⎧    ⎛               2  2                                                 
    ⎪    ⎜              z ⋅P                                                  
    ⎪    ⎜         ───────────────                                            
    ⎪    ⎜              ⎛  2     ⎞                                            
    ⎪    ⎜            4 ⎜ P     1⎟                                            
    ⎪    ⎜         4⋅Q ⋅⎜──── + ─⎟                                            
    ⎪    ⎜              ⎜   2   2⎟                                           2
    ⎪    ⎜              ⎝2⋅Q     ⎠    ⎛         z⋅P         ⎞               z 
    ⎪    ⎜  z⋅π⋅P⋅ℯ               ⋅erf⎜─────────────────────⎟          ───────
    ⎪    ⎜                            ⎜           __________⎟               ⎛ 
    ⎪    ⎜                            ⎜          ╱   2      ⎟             4 ⎜ 
    ⎪    ⎜                            ⎜   2     ╱   P     1 ⎟          4⋅Q ⋅⎜─
    ⎪    ⎜                            ⎜2⋅Q ⋅   ╱   ──── + ─ ⎟               ⎜ 
    ⎪    ⎜                            ⎜       ╱       2   2 ⎟               ⎝2
    ⎪    ⎜                            ⎝     ╲╱     2⋅Q      ⎠   z⋅π⋅P⋅ℯ       
    ⎪  Q⋅⎜- ───────────────────────────────────────────────── - ──────────────
    ⎪    ⎜                           __________                            ___
    ⎪    ⎜                          ╱   2                                 ╱   
    ⎪    ⎜                   2     ╱   P     1                     2     ╱   P
z ↦ ⎨    ⎜                2⋅Q ⋅   ╱   ──── + ─                  2⋅Q ⋅   ╱   ──
    ⎪    ⎜                       ╱       2   2                         ╱      
    ⎪    ⎝                     ╲╱     2⋅Q                            ╲╱     2⋅
    ⎪- ───────────────────────────────────────────────────────────────────────
    ⎪                                             3/2                         
    ⎪                                        2⋅z⋅π   ⋅P                       
    ⎪                                                                         
    ⎪                                                                         
    ⎪                                                                         
    ⎪                                                                         
    ⎪                                                                         
    ⎪                                                                         
    ⎪                                                                         
    ⎪                                                                         
    ⎪                                                                         
    ⎪                                                                         
    ⎪                                                                         
    ⎪                                                                         
    ⎪                                                                         
    ⎩                                                                         
... (426 lines in total)                                         

                         ⎞                                             ⎞
 ⎛      1          ⎞│   π⎟   │                 ⎛      1          ⎞│   π⎟
t⎜──────────────, ∞⎟│ < ─⎟ ∨ │periodic_argument⎜──────────────, ∞⎟│ < ─⎟
 ⎜          2      ⎟│   2⎟   │                 ⎜          2      ⎟│   2⎟
 ⎝polar_lift (P)   ⎠│    ⎠   │                 ⎝polar_lift (P)   ⎠│    ⎠

Почему SymPy дал такой странный результат и что я должен сделать, чтобы получить сжатую форму?

Answer 1

Объявите второй параметр (стандартное отклонение, «сигма») как положительный:

a = Normal('a', 0, Symbol('P', positive=True))
b = Normal('b', 0, Symbol('Q', positive=True))

Тогда результат, как и ожидалось:

               2     
             -z      
         ─────────── 
            2      2 
         2⋅P  + 2⋅Q  
     √2⋅ℯ            
z ↦ ─────────────────
            _________
           ╱  2    2 
    2⋅√π⋅╲╱  P  + Q  

Под капотом SymPy вычисляет неправильный интеграл, используя функцию Мейера G, которая включает в себя поднятие вычисления до некоторой римановой поверхности, а наличие отрицательных чисел для P или Q может привести к другой ветви; отсюда и сложный ответ. Математически, это не должно иметь большого значения для вычисления, так как сигма возводится в квадрат в показателе степени; но это имеет значение для успеха алгоритма.

Все это было бы спорным, если бы SymPy просто знал, как добавляются независимые нормали, но это не так; все такие вычисления выполняются путем прямой интеграции, часто выходящей за пределы реализованных методов интеграции.

Примечания:

1. Вы использовали true (объект правды SymPy), но Python True ожидается методом создания символа.
2. Лучше избегать объединения from sympy import * и from sympy.stats import *, поскольку это приводит к конфликтам имен: E - это число Эйлера в SymPy и - обозначение ожидаемого значения в sympy.stats.