There is another bug here with the order of `param_names`. This function should be:

    def eval(self, x, y, p):

        a = 0.5 * ((np.cos(p[5]) / p[3]) ** 2 +
                   (np.sin(p[5]) / p[4]) ** 2)
        b = 0.5 * (np.cos(p[5]) * np.sin(p[5]) *
                   (1. / p[3] ** 2 - 1. / p[4] ** 2))
        c = 0.5 * ((np.sin(p[5]) / p[3]) ** 2 +
                   (np.cos(p[5]) / p[4]) ** 2)

        return p[0] * np.exp(-(a * (x - p[1]) ** 2 +
                               b * (x - p[1]) * (y - p[2]) +
                               c * (y - p[2]) ** 2))

The eval was getting stddev's and mean's wrong.