Plotting confidence interval from P values (from DistributionFitTest) for 2 parameters

Question

I have a data set I'd like to adjust with 2 parameters, then compare that to another set of data. By using NMaximixe with DistributionFitTest, I'm able to determine the best fit of parameters. Great, now I'd like to do something like plot a ContourPlot of the 2D parameter space to show a 1-sigma, 2-sigma, and 3-sigma confidence interval for possible choices of parameters.

However, DistributionFitTest only returns P-values and test statistics. I am still a complete noob when it comes to stats, so I have no idea how to determine confidence intervals from this...

A quick example that you can copy and paste:

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];

(* 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];

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

(* 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"}]]

(* Now show a plot of the parameter space, where the peak should be around x = 1, y = 0.9 *)
ContourPlot[DistributionFitTest[(data1 - shift) scale, skd2, "PValue"], {shift, 0.8, 1.2}, {scale, 0.8, 1}]

I could spit out a test statistic for a grid of points instead (and try to fit lines to this), but I don't know how to calculate the confidence interval from P-values or test statistics. Is there a built-in method, or a way I could even brute force it?

Thanks in advance!

Answer 1

Here is a bootstrap approach where you don't have to assume any particular form of distribution. (And, yes, the code is made to show things step-by-step rather than being efficient. However, the 1,000 bootstrap samples takes only the blink of an eye.)

(* First, get an estimate of the shift and scale parameters *)
momScale = StandardDeviation[data2]/StandardDeviation[data1]
momShift = Mean[data1] - Mean[data2]/momScale
(* Or use the following *)
(* miScale=InterquartileRange[data2]/InterquartileRange[data1]
miShift=Median[data1]-Median[data2]/miScale *)

(* Now run 1,000 bootstrap samples *)
nBoot = 1000;
scale1 = 0;
scale2 = 0;
shift1 = 0;
shift2 = 0;
scaleXshift = 0;
Do[
 (* Take a bootstrap sample from each dataset *)
 d1 = RandomChoice[data1, NumPerSet];
 d2 = RandomChoice[data2, NumPerSet];

 (* Estimate parameters from the bootstrap sample *)
 bootScale = StandardDeviation[d1]/StandardDeviation[d2];
 bootShift = Mean[d1] - Mean[d2]/bootScale;

 (* Bump all of the counters used to estimate the standard errors *)
 scale1 = scale1 + bootScale;
 scale2 = scale2 + bootScale^2;
 shift1 = shift1 + bootShift;
 shift2 = shift2 + bootShift^2;
 scaleXshift = scaleXshift + bootScale*bootShift,
{i, nBoot}]

(* Standard error for scale *)
seScale = ((scale2 - scale1^2/nBoot)/(nBoot - 1))^0.5
(* Standard error for shift*)
seShift = ((shift2 - shift1^2/nBoot)/(nBoot - 1))^0.5
(* Correlation between estimators of shift and scale *)
corrScaleShift = ((scaleXshift - scale1 shift1/nBoot)/(nBoot - 1))/(seScale seShift)