code: def levinson_durbin(r, order):
      
      a = np.zeros(order + 1)
      e = np.zeros(order + 1)
  
      # k = 1
      a[0] = 1.0
      a[1] = - r[1] / r[0]
      e[1] = r[0] + r[1] * a[1]
      lam = - r[1] / r[0]
  
      for k in xrange(1, order):
          # update lambda
          lam = 0.0
          for j in xrange(k + 1):
              lam -= a[j] * r[k + 1 - j]
          lam /= e[k]
  
          # update a with U and V
          U = [1]
          U.extend([a[i] for i in xrange(1, k + 1)])
          U.append(0)
          V = [0]
          V.extend([a[i] for i in xrange(k, 0, -1)])
          V.append(1)
          a = np.array(U) + lam * np.array(V)
  
          # update e
          e[k + 1] = e[k] * (1.0 - lam * lam)
  
      return a, e[-1]
docs: from http://aidiary.hatenablog.com/entry/20120415/1334458954