Return to Revisions

With the help of comments on this and other stackexchange pages I managed to solve the problem of how to use custom distributions in things like CopulaDistribution (and other functions like RandomVariate, Expectation, etc), and given that it took me a couple of days’ hard slog I thought I’d share my discoveries with this community. Please excuse the flippant nature of the following – it started life as an email to a friend and colleague of mine… :)

(Also, please note that Random``DistributionVector, mentioned by Sasha in one of the linked posts above, appears to have been changed to Random\``Private``DistributionVector` in the production release of v9.0.1)

I needed a Mathematica version of Tukey’s “g-and-h” distribution, which is a transform of Z ~ N[0,1] ...

Xgh = A + B/g(Exp[g*Z] - 1)Exp[(h*Z^2)/2]

OK, you’d think just use

ghd = TransformedDistribution[Xgh, Distributed[Z, NormalDistribution[0,1]]

would you not?

Turns out that, even though the above is valid, it takes a very (VERY) long time to compute CDF[ghd, X], or even Quantile[ghd, p]. I tried a number of tricks including Compile, taking a Taylor series, etc, but I’ve concluded (reluctantly) it’s not viable to use TransformedDistribution for g-and-h.

OK then, you’d think just use the g-and-h cdf prob in

ghd = ProbabilityDistribution[{“CDF”, prob}, {x, -Infinity, Infinity}] 

would you not?

Turns out that a symbolic form of the pdf and/or cdf is problematic, and using

prob = CDF[NormalDistribution[0, 1], (Z /. FindRoot[Xgh == x, {Z,0}])] 

doesn’t cut it with Mathematica. It wants something symbolic inside ProbabilityDistribution.

So, I reverted back to my “hand-rolled” (custom) distribution that I’d written before the new fancy enhanced Distributions came out in v8 & v9. It already had all the std Mathematica things like CDF (implemented as the FindRoot of Xgh), PDF, InverseCDF, Random, Quantile, Mean, Variance, etc, via the TagSet feature (ie, defining “up values” of the distribution name). Eg...

GandHDistribution/: InverseCDF[GandHDistribution[A_, B_, g_, h_], fraction_] :=
  Xgh[A, B, g, h, Quantile[NormalDistribution[0, 1], fraction]] /; 
    DistributionParameterQ[GandHDistribution[A, B, g, h]]

But it turns out that there are a lot of undocumented test functions used by CopulaDistribution that you also need to define. Also, some of the pre-v8 internals had changed, e.g., ParameterQ is now DistributionParameterQ. I found out what was needed by Unprotect[CopulaDistribution] and then removing the ReadProtected attribute so that I can see all the required up-value definitions. They appear to be (aside from the “well-known” things like InverseCDF and DistributionParameterQ mentioned above)...

DistributionDomain

DistributionParameterAssumptions

GandHDistribution/: 
  Random`Private`
    DistributionVector[GandHDistribution[A_, B_, g_, h_],
                       n_Integer, prec_?Positive] :=
      Xgh[A, B, g, h, RandomVariate[NormalDistribution[0, 1], n, 
        WorkingPrecision -> prec]];

GandHDistribution/: 
 Statistics`CopulaDistributionDump` 
   UnivariateDistributionListQ[GandHDistribution[A_, B_, g_, h_]] := True;

GandHDistribution/: 
  Statistics`Library`
    ContinuousUnivariateDistributionQ[GandHDistribution[A_, B_, g_, h_]] := True;

GandHDistribution/: 
  Statistics`Library`
    DiscreteUnivariateDistributionQ[GandHDistribution[A_, B_, g_, h_]] := False;

GandHDistribution/: 
  Statistics`Library`
    ContinuousMultivariateDistributionQ[GandHDistribution[A_, B_, g_, h_]] := False;

GandHDistribution/: 
  Statistics`Library`
    DiscreteMultivariateDistributionQ[GandHDistribution[A_, B_, g_, h_]] := False;

GandHDistribution/: 
  Statistics`Library`
    DistributionNParameterQ[GandHDistribution[A_, B_, g_, h_]]:=
     DistributionParameterQ[GandHDistribution[A, B, g, h]];

Now it works like a dream – the CopulaDistribution function is fantastic.