1

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!

Matt Stein
  • 145
  • 8
  • By now I think you know what my comment will be: Get thee to a Statistician! Your approach is imaginative but certainly not commonly accepted statistical practice. – JimB Oct 06 '17 at 18:23

1 Answers1

2

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)
JimB
  • 41,653
  • 3
  • 48
  • 106
  • Awesome, that's quick and easy. Now... plotting this... In theory, I should be able to position the "peak" at (momShift,momScale) on a 2D plot of the parameter space, and see rings go out at 1/2/3-sigma, but the exact shape/rotation of such a contour plot eludes me. For example, is it rotated at all? – Matt Stein Oct 06 '17 at 20:28
  • Forget 1/2/3-sigma: at best that applies to the one dimensional case. One is assuming/hoping that a bivariate normal distribution will closely approximate the bivariate distribution of the estimators of the scale and shift parameters. Search this site for "bivariate normal contours" or see https://mathematica.stackexchange.com/questions/153511/1-sigma-contour-plot-for-a-two-dimensional-gaussian/153516#153516. – JimB Oct 06 '17 at 20:46
  • If you can bear with me a while longer: Forget the p-value approach, too. Your problem as stated involves "estimation" rather than "hypothesis testing". – JimB Oct 06 '17 at 20:49
  • I am 100% OK abandoning the P-value approach, and you're right, I'm just trying to get good estimations. I will refer to your other link for more help with the contours. Thank you again!! – Matt Stein Oct 06 '17 at 20:53