2016-04-16 2 views
1

Для проекта я хочу использовать sympy для построения и вычисления максимальной вероятности гауссовского распределения для дискретного числа точек данных. Метод, который я следую, находится по адресу mathworld.Как использовать массив python в символическом выражении?

Но я столкнулся с проблемой, когда попытался использовать массив в символическом выражении с Product и/или Sum. Ниже приведен упрощенный вариант моих предыдущих попыток.

В записной книжке в Yupyter Anaconda, я создать массив питона, скажем x:

N = 10 
x = range(N) 

И я хочу использовать x в символическое выражение в sympy следующим образом:

from sympy import *, Symbol 
i = Symbol('i', integer=True) 
mu = Symbol('mu') 

s = Sum((x[i]-mu)**2, (i, 0, N-1)) 

Но это не работает, поскольку оценка ячейки приводит к:

TypeError       Traceback (most recent call last) 
<ipython-input-1-19c174235872> in <module>() 
     7 mu = Symbol('mu') 
     8 
----> 9 s = Sum((x[i]-mu)**2, (i,0,N-1)) 

TypeError: list indices must be integers, not Symbol 

Еще одна попытка:

X = MatrixSymbol(X, 1, N) # No clue how to convince `sympy` to use 1-dimensional arrays using only one index. 

s = Sum((X[0,i]-mu)**2, (i,0,N-1)) 
s.doit() 

Дает:

(-mu + X[0, 0])**2 + (-mu + X[0, 1])**2 + (-mu + X[0, 2])**2 + (-mu + X[0, 3])**2 + (-mu + X[0, 4])**2 + (-mu + X[0, 5])**2 + (-mu + X[0, 6])**2 + (-mu + X[0, 7])**2 + (-mu + X[0, 8])**2 + (-mu + X[0, 9])**2 

вид работ, но как получить реальные значения x в этом символическое выражение, т.е. заменяя каждый из этих X[0,i] для значения от x[i]?

Еще одна попытка:

X = Matrix(1,N, range(N)) 

s = Sum((X[i]-mu)**2, (i, 0, N-1)) 
s.doit() 

Сейчас питон/SymPy очень несчастна:

--------------------------------------------------------------------------- 
IndexError        Traceback (most recent call last) 
<ipython-input-7-6d106fb975e1> in <module>() 
     1 X = Matrix(1,N, range(N)) 
     2 
----> 3 s = Sum((X[i]-mu)**2, (i, 0, N-1)) 
     4 s.doit() 

/Users/twan/anaconda/lib/python2.7/site-packages/sympy/matrices/dense.pyc in __getitem__(self, key) 
    94    if isinstance(key, slice): 
    95     return self._mat[key] 
---> 96    return self._mat[a2idx(key)] 
    97 
    98  def __setitem__(self, key, value): 

/Users/twan/anaconda/lib/python2.7/site-packages/sympy/matrices/matrices.pyc in a2idx(j, n) 
    4412    j = j.__index__() 
    4413   except AttributeError: 
-> 4414    raise IndexError("Invalid index a[%r]" % (j,)) 
    4415  if n is not None: 
    4416   if j < 0: 

IndexError: Invalid index a[i] 

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

EDIT:

Я должен отметить, что слагаемым должен был быть разработан в (E**(-((x-mu)**2)/(2 * s**2)))/(s * sqrt(2 * pi)). Возрождение делает его несколько иной проблемой.

Хотя решение @ unutbu не сработало для меня, пытаясь использовать это предложение, указал мне, что я должен сдерживать домены до действительных чисел.

@ предложение Зефир в сделал работу, полное решение теперь:

from sympy import symbols, E, pi, sqrt, init_printing 
from sympy import diff, IndexedBase 
from sympy.solvers import solve 

x, mu = symbols('x mu', real=True) 
sigma = symbols('sigma', real=True, positive=True) 

bell = (E**(-((x-mu)**2)/(2 * sigma**2)))/(sigma * sqrt(2 * pi)) 

def likelihood(factor, xs): 
    return np.prod([factor.subs(x,i) for i in xs]) 

def loglikelihood(factor, xs): 
    return expand_log(log(likelihood(factor, xs))) 

N = 3 
X = IndexedBase('X') 
Xs = [X[i] for i in range(N)] 

solve(diff(loglikelihood(gauss,Xs), mu).subs(sigma, 1), mu) 

Большое спасибо @Marshmellow и @unutbu.

ответ

1

Поскольку ваш список x является числовым, для его обработки вам не требуется символическая сумма.Просто просуммировать список (xi-mu)**2 с помощью языка Python список понимание:

from sympy import * 
N = 10 
x = range(N) 
mu = Symbol('mu') 
s = sum([(xi-mu)**2 for xi in x]) 
print(s) 
print(s.diff(mu))  # to show this is a symbolic expression 

Выход:

mu**2 + (-mu + 1)**2 + (-mu + 2)**2 + (-mu + 3)**2 + (-mu + 4)**2 + (-mu + 5)**2 + (-mu + 6)**2 + (-mu + 7)**2 + (-mu + 8)**2 + (-mu + 9)**2 
20*mu - 90 
+0

Спасибо @Marshmallow, ваше предложение использовать сумму или продукт над пониманием списка поставил меня на правильный путь. См. Мое редактирование. – nanitous

1

Вы можете использовать IndexedBase представлять массив, содержащий элементы.

X = sy.IndexedBase('X') 
s = sy.Sum((X[i]-mu)**2, (i, 0, N-1)) 

можно использовать lambdify заменить символы SymPy с NumPy массивы. В этом случае мы хотим заменить X[i] на значения из массива NumPy.

В настоящее время lambdifycan not be applied to IndexedBase объектов. Но он может применяться к DeferredVector с. Например:

import sympy as sy 
import numpy as np 

i = sy.Symbol('i', integer=True) 
mu = sy.Symbol('mu') 

N = 10 
X = sy.IndexedBase('X') 
s = sy.Sum(sy.exp((X[i]-mu)**2), (i, 0, N-1)) 

f = sy.lambdify(sy.DeferredVector('X'), s, 'sympy') 
x = np.arange(N) 
print(f(x)) 

печатает

exp(mu**2) + exp((-mu + 1)**2) + exp((-mu + 2)**2) + exp((-mu + 3)**2) + exp((-mu + 4)**2) + exp((-mu + 5)**2) + exp((-mu + 6)**2) + exp((-mu + 7)**2) + exp((-mu + 8)**2) + exp((-mu + 9)**2) 

Следует отметить, что, поскольку f(x) еще выражение SymPy, я использовал 'sympy' в качестве третьего аргумента lambdify так, что sy.exp не заменяется числовой exp функции.

+0

ваше решение, похоже, имеет проблемы, когда я использую 'sy.exp (X [i] -mu) ** 2)' как слагаемое. Sympy жалуется, что не может печатать поплавок! Это не слишком удивительно, поскольку я могу себе представить, что решение во всей его общности представлено в Sympy как выражение мнимого числа. – nanitous

+0

Я не уверен, как ограничить домен индексированной базы 'X' реальными числами, как вы можете, когда объявляете обычные символы с помощью символа. – nanitous

+0

Поскольку 'mu' по-прежнему является символом,' f (x) 'по-прежнему является символом. Поэтому вы не хотите, чтобы 'sy.exp' заменялся числовой функцией' exp'. Используйте 'f = sy.lambdify (sy.DeferredVector ('X'), s, 'sympy')', чтобы предотвратить замену 'sy.exp' на числовую функцию exp. – unutbu

Смежные вопросы