# Using DistributionFitTest on custom distributions in Mathematica 8

I originally posted this question on Stack Overflow, but I didn't get any answers, and I'm hoping to have better luck here.

I'm trying to compute the goodness-of-fit of a bi-modal Gaussian distribution. To do this, Mathematica seems to require a symbolic distribution function to which to compare. Because such a bi-modal distribution is not a stock distribution, I'm trying to define one.

The obvious use of

MixtureDistribution[{fs, (1-fs),
{NormalDistribution[μS, σS], NormalDistribution[μL, σL]}]


generates a distribution that can be plotted, but the analysis used by DistributionFitTest[] fails.

This topic has been addressed in previous questions in discussions between @Sasha and @Jagra:

http://stackoverflow.com/questions/6361647/distributionfittest-for-custom-distributions-in-mathematica

http://stackoverflow.com/questions/9202531/minimizing-nexpectation-for-a-custom-distribution-in-mathematica?lq=1

but I was unable to find a resolution that enabled the use of

DistributionFitTest[data,dist,"HypothesisTestData"]


when dist is not a built-in distribution type.

Because the distribution I'm modeling is composed of simple pieces, describing the properties of the distribution is not too difficult, and I have attempted to describe as many features as I know in order to create a well defined distribution that Mathematica 8 would recognize as one of its own. My attempt to define every parameter I can think of follows:

modelDist /:
PDF[modelDist[fS_, μS_, σS_, μL_, σL_], x_] :=
PDF[MixtureDistribution[{fS, 1 - fS}, {NormalDistribution[μS, σS], NormalDistribution[μL, σL]}], x];

modelDist /:
CDF[modelDist[fS_, μS_, σS_, μL_, σL_], x_] :=
CDF[MixtureDistribution[{fS, 1 - fS}, {NormalDistribution[μS, σS], NormalDistribution[μL, σL]}], x];

modelDist /:
DistributionDomain[modelDist[fS_, μS_, σS_, μL_, σL_]] :=
Interval[{-Infinity, Infinity}];

modelDist /:
RandomDistributionVector[modelDist[fS_, μS_, σS_, μL_, σL_], n_, prec_] :=
RandomVariate[MixtureDistribution[{fS, 1 - fS}, {NormalDistribution[μS, σS], NormalDistribution[μL, σL]}], n, WorkingPrecision -> prec];

modelDist /:
DistributionParameterQ[modelDist[fS_, μS_, σS_, μL_, σL_]] :=
!TrueQ[Not[Element[{fS, μS, σS, μL, σL}, Reals] && fS > 0 && fS < 1 && σS > 0 && σL > 0]];

modelDist /:
DistributionParameterAssumptions[modelDist[fS_, μS_, σS_, μL_, σL_]] :=
Element[{fS, μS, σS, μL, σL}, Reals] && fS > 0 && fS < 1 && σS > 0 && σL > 0;

modelDist /:
MomentGeneratingFunction[modelDist[fS_, μS_, σS_, μL_, σL_], t_] :=
fS E^(t μS + (t^2 σS^2)/2) + (1 - fS) E^(t μL + (t^2 σL^2)/2);

modelDist /:
CharacteristicFunction[modelDist[fS_, μS_, σS_, μL_, σL_], t_] :=
fS E^(I t μS + (t^2 σS^2)/2) + (1 - fS) E^(I t μL + (t^2 σL^2)/2)

modelDist /:
Moment[modelDist[fS_, μS_, σS_, μL_, σL_], n_] :=
Piecewise[{{fS*σS^n*(-1 + n)!!*Hypergeometric1F1[-(n/2), 1/2, -(μS^2/(2*σS^2))] + (1 - fS) * σL^n*(-1 + n)!! * Hypergeometric1F1[-(n/2), 1/2, -(μL^2/(2*σL^2))], Mod[n, 2] == 0}}, μS*σS^(-1 + n)*n!!* Hypergeometric1F1[(1 - n)/2, 3/2, -(μS^2/(2*σS^2))] + (1 - fS) * μL*σL^(-1 + n)*n!! * Hypergeometric1F1[(1 - n)/2, 3/2, -(μL^2/(2*σL^2))]];

modelDist /:
Mean[modelDist[fS_, μS_, σS_, μL_, σL_]] :=
fS μS + (1 - fS) μL

modelDist /:
Expectation[expr_, x_ \[Distributed] modelDist[fS_, μS_, σS_, μL_, σL_]] :=
fS*Expectation[expr, x \[Distributed] NormalDistribution[μS, σS]] + (1 - fS)*Expectation[expr, x \[Distributed] NormalDistribution[μL, σL]]


Everything seems to work up through the definition of Expectation, which throws

TagSetDelayed::tagpos: Tag modelDist in Expectation[expr_, x_\[Distributed]modelDist[fS_, \μS_, \σS_, \μL_, \σL_]] is too deep for an assigned rule to be found.

I don't know that having a definition for the expectation will magically make everything work, but it's the next step to to try, as having the expectation allows computation of the variance, and for all I know, that is the last tag that I need to define. Is there a syntax that will properly define this Expectation in such a way that the expression will pass straight from my modelDist to its constituent NormalDistributions?

-
I wonder if you could run the distribution fit tests individually? It might help narrow down what doesn't work. –  Jagra Dec 31 '12 at 16:39

Not really an answer but too long for a comment. This works for me in versions 8.0.4 and 9.0.0.

dist =
modelDist[fS, μS, σS, μL, σL] /.
FindInstance[
DistributionParameterAssumptions[
modelDist[fS, μS, σS, μL, σL]],
{fS, μS, σS, μL, σL}, Reals][[1]];

DistributionFitTest[RandomVariate[dist, 100], dist]

(* 0.846971 *)


It isn't necessary for Expectation to work for DistributionFitTest to do its job. DistributionFitTest does need CDF to be defined which you have done.

-
This generates a specific instance of modelDist, but it doesn't find the parameters fS,\[Mu]S, etc., that best fit a data sample. This is possible using htd=DistributionFitTest[RandomVariate[dist,100],NormalDistribution[\[Mu],\[Sigm‌​a]], "HypothesisTestData", even if the uni-modal fit doesn't match the bimodal distribution. htd will contain best fit values for \[Sigma] and \[Mu]. However, I can't figure out how to plug modelDist directly into DistributionFitTest so that I can get estimates and confidence intervals for fS,\[Mu]S,\[Sigma]S,\[Mu]L, and \[Sigma]L. –  KDN Dec 31 '12 at 14:58
This is much more complicated. When estimating the parameters from the data itself the test is no longer valid (unless a correction factor has been added for that distribution with estimated parameters which it hasn't). You can still use EstimatedDistribution to get the maximum likelihood estimate if that is what you are after but the test on the resulting distribution is still not valid. You will need to resort to other methods such as quantile plots etc. –  Andy Ross Dec 31 '12 at 16:58
What I'm after in this instance is the "plug-n-play"-ness of being able to drop a set of data into DistributionFitTest against a model distribution with variable parameters. If the distribution I was trying to fit was more complicated, I might give up, but it's just the sum of two gaussians. It seems this should be very straightforward, and I'd happily specify all behaviors by hand in order to get DistributionFitTest` to work in this fashion, as it does with the "native" Mathematica distributions. However, I don't know what is needed in this regard. –  KDN Dec 31 '12 at 18:39