FindMaximum for LogLikelihood fails for multiple parameters

Question

I need to fit a non-analytic data set to another data set using at least 2 parameters: shift and scale. For my 2 data sets (data1 and data2), I will to apply these 2 parameters to data1 such that it matches data2. This works perfectly for only one of those parameters, but when I try to do both at the same time, it fails. The 2nd parameter is never actually fit...

I will illustrate this by creating a "double hump" distribution using two normal distributions.

Clear["Global`*"]

NumPerSet = 1000;

(* This is the first "double-hump" distribution *)
set1 = RandomVariate[NormalDistribution[6, 1.2], NumPerSet];
set2 = RandomVariate[NormalDistribution[11, 1.2], NumPerSet];
data1 = Join[set1, set2];

Then I do the same thing again, but with the peaks shifted to the left by 1:

(* This is like data1 but "shifted" to the left by 1 *)
set3 = RandomVariate[NormalDistribution[5, 1.2], NumPerSet];
set4 = RandomVariate[NormalDistribution[10, 1.2], NumPerSet];
data2 = Join[set3, set4];

Then, I apply the scale factor, where each data point in the set has been reduced by 10%:

(* Now apply a scale factor to the data2 set *)
data2 = data2 0.9;

To visualize this, I can plot:

(* Create some plottable distributions *)
skd1 = KernelMixtureDistribution[data1];
skd2 = KernelMixtureDistribution[data2];

(* Show each distribution and the associated histogram *)
Show[Plot[{PDF[skd1, x], PDF[skd2, x]}, {x, 0, 15}, Frame -> True, PlotLabel -> "The Setup"], Histogram[{data1, data2}, Automatic, "PDF", ChartLegends -> {"data1", "data2"}]]

Which returns:

Great, now to solve for the parameters:

(*Mean of distribution*)
\[Mu] = Expectation[x, x \[Distributed] skd2];

(*Estimate a by method of moments*)
mom = Abs[Mean[data1] - \[Mu]];

(*Estimate a by maximum likelihood*)
logL = LogLikelihood[skd2, (data1 - shift) scale];
fmout = FindMaximum[logL, {{shift, mom}, {scale, .9}}, AccuracyGoal -> 20]

Print["fmout = ", fmout]
Print["Max Likelihood value for:"]
Print["  Shift = ", Part[Part[Part[fmout, 2], 1], 2]]
Print["  Scale = ", Part[Part[Part[fmout, 2], 2], 2]]
Print["  MLE   = ", Part[fmout, 1]]

(* Now show the actual adjustment on data1 which supposedly fits data2 the best *)
Show[Plot[{PDF[skd2, x]}, {x, 0, 15}, Frame -> True, PlotLabel -> "Fitted Values"], Histogram[{data2, (data1 Part[Part[Part[fmout, 2], 2], 2]) -     Part[Part[Part[fmout, 2], 1], 2]}, Automatic, "PDF", ChartLegends -> {"data1", "data2"}]]

(* Now show what we already know is the true shift and scale *)
Show[Plot[{PDF[skd2, x]}, {x, 0, 15}, Frame -> True, PlotLabel -> "Already Known Values"], Histogram[{data2, (data1 - 1) .9}, Automatic, "PDF"]]

What consistently comes back is that the first parameter (shift) is fit, but the second parameter is just stuck at whatever I put in for the starting value (in this case, 0.9).

For example, in one case, the shift was 1.76895, scale was 0.9, and the plot:

And by applying the "already known" values of shift = 1, scale = 0.9, I get:

Obviously a much better fit... Can anyone help me understand why this is failing?

Also, if I try to plot the MLE via ContourPlot, I get something that looks like the whole thing is just borked (shift on X-axis, scale on Y-axis, red dot indicates the "fitted" parameters):

Thanks!

edit: Code for the contour plot:

Show[ContourPlot[logL, {shift, 0, 2}, {scale, 0.8, 1.0}], ListPlot[{{Part[Part[Part[fmout, 2], 1], 2], Part[Part[Part[fmout, 2], 2], 2]}}, PlotStyle -> Red]]

Answer 1

If you really know that your distributions differ only by a shift and scale, then either the method of moments or the use of median and interquartile range will work much better and faster.

Method of moments:

momScale = StandardDeviation[data2]/StandardDeviation[data1]
momShift = Mean[data1] - Mean[data2]/momScale

Median and interquartile range:

mirScale = InterquartileRange[data2]/InterquartileRange[data1]
mirShift = Median[data1] - Median[data2]/mirScale

And by "better" I mean that the standard errors are much smaller than the two approaches you've described above. (And you don't have to put in possibly unrealistic bounds on the parameters.) Of course, I'm assuming that for the method of moments the first two moments exist for the underlying distribution. Otherwise, I'm not assuming any particular shape of distribution.