ENH: show the fit on top of data points.

As discussed on the mailinglist the default for ModelResult.plot_fit()
has been changed to show the fit on top of the data points.

(see https://groups.google.com/forum/#!topic/lmfit-py/rekwbFVJy2c)

lmfit/model.py

@@ -1649,11 +1649,6 @@ def plot_fit(self, ax=None, datafmt='o', fitfmt='-', initfmt='--',
                                 **{independent_var: x_array_dense}),
                 initfmt, label='init', **init_kws)
 
-        ax.plot(
-            x_array_dense,
-            self.model.eval(self.params, **{independent_var: x_array_dense}),
-            fitfmt, label='best-fit', **fit_kws)
-
         if yerr is None and self.weights is not None:
             yerr = 1.0/self.weights
         if yerr is not None:
@@ -1662,6 +1657,11 @@ def plot_fit(self, ax=None, datafmt='o', fitfmt='-', initfmt='--',
         else:
             ax.plot(x_array, self.data, datafmt, label='data', **data_kws)
 
+        ax.plot(
+            x_array_dense,
+            self.model.eval(self.params, **{independent_var: x_array_dense}),
+            fitfmt, label='best-fit', **fit_kws)
+
         ax.set_title(self.model.name)
         if xlabel is None:
             ax.set_xlabel(independent_var)