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 MMA version of Tukey’s “g-and-h” distribution, which is a transform of Z ~ N[0,1] ...
Xgh = A + B/g(Exp[gZ] - 1)Exp[(hZ^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, eg. 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 = Z /. FindRoot[Xgh == x, {Z,0} ] doesn’t cut it with MMA, 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 MMA things like CDF, PDF, InverseCDF(implemented as the FindRoot[ ] of Xgh), 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, eg. 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 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.
Cheers, Doug.