def N(o_t, mu_j, sigma_j):
    p = len(o_t)
    o_t = matrix(o_t).T # px1 vector
    mu_j = matrix(mu_j).T # px1 vector
    sigma_j = matrix(sigma_j) # pxp matrix
    det = linalg.det(sigma_j)
    normalizer = 1.0 / (pow(2*pi,p/2)*sqrt(det))
    a = o_t - mu_j
    d = normalizer * exp(-0.5 * a.T * sigma_j.I * a)
    return d