4

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:

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 /:
Random`DistributionVector[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 `NormalDistribution`

s?

(And if this entirely the wrong way to go about this, some advice to that effect would be appreciated.)

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