How do I use a python array in a symbolic expression?

For a project I want to use sympy to construct and compute the maximum likelihood of a gaussian distribution for a discrete number of data points. The method I follow can be found at mathworld.

But I ran into trouble when I attempted to use an array in a symbolic expression with Product and/or Sum. Below is a simplified version of my earlier attempts.

In notebook in Anaconda's Yupyter, I create a python array, say x:

N = 10
x = range(N)

And I want to use x in a symbolic expression in sympy as follows:

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

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

But this doesn't work, because evaluating the cell results in:

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

Another attempt:

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()

Gives:

(-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

kind of works, but how to get the real values for x into that symbolic expression, i.e. substituting each of these X[0,i] for the value of x[i]?

Another attempt:

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

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

Now is python/sympy very unhappy:

---------------------------------------------------------------------------
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                 retu self._mat[key]
---> 96             retu 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]

I've no clue what to try else and I'm stuck with sympy here. I'm curious whether I ran into a limitation of sympy for an very important computation in statistics.