(*****************************************************************************
Conditions of Use

Of the Mathematica package JohnsonDistributions.m

by Sergej V. Aksenov (http://aksenov.freeshell.org/)


This software (file JohnsonDistributions.m) is copyright by Sergej V. Aksenov 
and Michael A. Savageau.

This software and its constituent parts come with no warranty, whether expressed 
or implied, that it is free of errors or suitable for any specific purpose.  It 
must not be used to solve any problem, whose incorrect solution could result in 
injury to a person, institution, or property.  It is user's own risk to use this 
software or its parts and the authors disclaim all liability for such use.

This software is distributed ``as is''. In particular, no maintenance, support, 
troubleshooting, or subsequent upgrade is implied.

The use of this software must be acknowledged in any publication that contains 
results obtained with it; please contact the author for the citation.

The free use of this software and its parts is restricted for research 
purposes; commercial use require permission and licensing from Sergej V. 
Aksenov and Michael A. Savageau.
******************************************************************************)



(*:Mathematica Version: 4.1 *)

(*:Package Version: 1.21 *)

(*:Name: JohnsonDistributions` *)

(*:Context: JohnsonDistributions`JohnsonDistributions` *)

(*:Title: Bounded and unbounded Johnson distributions *)

(*:Author: Sergej V. Aksenov (http://aksenov.freeshell.org/) *)

(*:History:
   Created: August, 2001.
   Last modified: April 24, 2002: added some moments of Johnson distributions.
   Last modified: April 30, 2002: added ParameterQ and passing options to NIntegrate. 
   Last modified May 1, 2002: added exact calculation of unbounded distribution moments.
   v1.21 Last modified 9 May, 2002: implemented proper FilterOptions with NIntegrate.
*)

(*:Copyright: Copyright 2001-2, Sergej V. Aksenov and Michael A. Savageau. *)

(*:Reference: Usage messages. *)

(*:Summary:
This package provides properties and functionals of continuous bounded (S_B) and unbounded (S_U) Johnson distributions. *)

(*:Keywords: 
Continuous distribution, Johnson distributions. *)

(*:Limitations: 
CharacteristicFunction, ExpectedValue are not implemented. 
*)

(*:Note: 
Meanm Variance, StandardDeviation, Skewness, Kurtosis for unbounded distribution are calculated using exact formulas.
Mean, Variance, StandardDeviation, Skewness, Kurtosis for bounded distribution are calculated via NIntegrate.  To calculate moments and their functions to higher precision one may need to use SetOptions[NIntegrate, PrecisionGoal->prec], or pass this option to a statistical function directly, it will be transmitted to NIntegrate through FilterOptions.
*)

(*:Sources:
  N. L. Johnson, Systems of frequency curves generated by methods of translation, Biometrika 36 (1-2): 149-176, 1949.
*)



BeginPackage["JohnsonDistributions`JohnsonDistributions`","Statistics`Common`DistributionsCommon`","Statistics`DescriptiveStatistics`","Statistics`NormalDistribution`","Utilities`FilterOptions`"]



Unprotect[UnboundedJohnsonDistribution, BoundedJohnsonDistribution]



(* Distributions. *)

(* If z=g+d*ArcSinh[x] has standard Normal distribution, then x is unbounded Johnson (S_U) random variable. *)

UnboundedJohnsonDistribution::usage="UnboundedJohnsonDistribution[g,d] represents unbounded Johnson distribution with shape parameters g and d."

(* If z=g+d*Log[x/(1-x)] has standard Normal distribution, then x is bounded Johnson random variable. *)

BoundedJohnsonDistribution::usage="BoundedJohnsonDistribution[g,d] represents bounded Johnson distribution with shape parameters g and d."



Begin["`Private`"]



(* Functionals of Unbounded distribution. *)

UnboundedJohnsonDistribution/: ParameterQ[UnboundedJohnsonDistribution[g_, d_]] := And[ If[FreeQ[N[p], Complex], True, Message[UnboundedJohnsoDistribution::realparm, g]; False], If[FreeQ[N[q], Complex], True, Message[UnboundedJohnsoDistribution::realparm, d]; False],  If[N[d] > 0, True, Message[UnboundedJohnsonDistribution::posparm, d]; False, True]]

UnboundedJohnsonDistribution::realparm = "Parameter `1` is expected to be real."
UnboundedJohnsonDistribution::posparm = "Parameter `1` is expected to be positive."

UnboundedJohnsonDistribution/: Domain[UnboundedJohnsonDistribution[g_, d_]] := Interval[{-Infinity, Infinity}]

UnboundedJohnsonDistribution/: PDF[UnboundedJohnsonDistribution[g_, d_], x_] := 1./Sqrt[2 Pi] Exp[-(g+d ArcSinh[x])^2/2] d/Sqrt[1+x^2] /; ParameterQ[UnboundedJohnsonDistribution[g,d]]

UnboundedJohnsonDistribution/: CDF[UnboundedJohnsonDistribution[g_, d_], x_] := 0.5 (Erf[(g+d ArcSinh[x])/Sqrt[2.]] + 1.0)  /; ParameterQ[UnboundedJohnsonDistribution[g,d]]

UnboundedJohnsonDistribution/: Quantile[UnboundedJohnsonDistribution[g_, d_], fraction_] := Sinh[(Quantile[NormalDistribution[0, 1], fraction] - g)/d] /; QuantileQ[fraction] && ParameterQ[UnboundedJohnsonDistribution[g,d]]

UnboundedJohnsonDistribution/: Random[UnboundedJohnsonDistribution[g_, d_]] := Quantile[UnboundedJohnsonDistribution[g, d], Random[]]

UnboundedJohnsonDistribution/: RandomArray[UnboundedJohnsonDistribution[g_, d_], dim_] := Table[Quantile[UnboundedJohnsonDistribution[g, d], Random[]], {dim}] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]

(* UnboundedJohnsonDistribution/: Mean[UnboundedJohnsonDistribution[g_, d_], opts___?OptionQ] := Module[{x, ret}, With[{optn = FilterOptions[NIntegrate, Sequence@@Flatten[{opts}]]}, ret = NIntegrate[ x PDF[UnboundedJohnsonDistribution[g, d], x], {x, -Infinity, Infinity}, optn]]; Return[ret]]  /; ParameterQ[UnboundedJohnsonDistribution[g,d]]
*)

(* Direct calculation of moments. *)

UnboundedJohnsonDistribution/: Mean[UnboundedJohnsonDistribution[g_, d_], opts___?OptionQ] := Module[{om}, om=Exp[1/d^2]; -Sqrt[om] Sinh[g/d] //N]  /; ParameterQ[UnboundedJohnsonDistribution[g,d]]

(* UnboundedJohnsonDistribution/: Variance[UnboundedJohnsonDistribution[g_, d_], opts___?OptionQ] := Module[{mu,x}, mu = Mean[UnboundedJohnsonDistribution[g, d], opts]; NIntegrate[ (x-mu)^2 PDF[UnboundedJohnsonDistribution[g, d], x], {x, -Infinity, Infinity}, opts]]  /; ParameterQ[UnboundedJohnsonDistribution[g,d]]
*)

UnboundedJohnsonDistribution/: Variance[UnboundedJohnsonDistribution[g_, d_], opts___?OptionQ] :=  Module[{om}, om=Exp[1/d^2]; (om-1)/2 (om Cosh[2 g/d]+1) //N] /; ParameterQ[UnboundedJohnsonDistribution[g,d]]

(* UnboundedJohnsonDistribution/: StandardDeviation[UnboundedJohnsonDistribution[g_, d_], opts___?OptionQ] := Module[{mu,x}, mu = Mean[UnboundedJohnsonDistribution[g, d], opts]; Sqrt[ NIntegrate[ (x-mu)^2 PDF[UnboundedJohnsonDistribution[g, d], x], {x, -Infinity, Infinity}, opts]]] /; ParameterQ[UnboundedJohnsonDistribution[g,d]]
*)

UnboundedJohnsonDistribution/: StandardDeviation[UnboundedJohnsonDistribution[g_, d_], opts___?OptionQ] :=  Module[{om}, om=Exp[1/d^2]; Sqrt[(om-1)/2 (om Cosh[2 g/d]+1)] //N] /; ParameterQ[UnboundedJohnsonDistribution[g,d]]

(* UnboundedJohnsonDistribution/: Skewness[UnboundedJohnsonDistribution[g_, d_], opts___?OptionQ] := Module[{sig,mu,x}, mu = Mean[UnboundedJohnsonDistribution[g, d], opts]; sig = Sqrt[NIntegrate[ (x-mu)^2 PDF[UnboundedJohnsonDistribution[g, d], x], {x, -Infinity, Infinity}, opts]]; 1/sig^3 NIntegrate[ (x-mu)^3 PDF[UnboundedJohnsonDistribution[g, d], x], {x, -Infinity, Infinity}, opts]] /; ParameterQ[UnboundedJohnsonDistribution[g,d]]
*)

UnboundedJohnsonDistribution/: Skewness[UnboundedJohnsonDistribution[g_, d_], opts___?OptionQ] := Module[{sig,mu3,om,ombig}, om=Exp[1/d^2]; ombig=g/d; mu3 = -Sqrt[om]/4 (om-1)^2 (om (om+2) Sinh[3 ombig]+3 Sinh[ombig]); sig = Sqrt[(om-1)/2 (om Cosh[2 ombig]+1)]; mu3/sig^3 //N] /; ParameterQ[UnboundedJohnsonDistribution[g,d]]

(* UnboundedJohnsonDistribution/: Kurtosis[UnboundedJohnsonDistribution[g_, d_], opts___?OptionQ] := Module[{sig,mu,x}, mu = Mean[UnboundedJohnsonDistribution[g, d], opts]; sig = Sqrt[NIntegrate[ (x-mu)^2 PDF[UnboundedJohnsonDistribution[g, d], x], {x, -Infinity, Infinity}, opts]]; 1/sig^4 NIntegrate[ (x-mu)^4 PDF[UnboundedJohnsonDistribution[g, d], x], {x, -Infinity, Infinity}, opts]] /; ParameterQ[UnboundedJohnsonDistribution[g,d]]
*)

UnboundedJohnsonDistribution/: Kurtosis[UnboundedJohnsonDistribution[g_, d_], opts___?OptionQ] := Module[{sig,mu4,om,ombig}, om=Exp[1/d^2]; ombig=g/d; sig = Sqrt[(om-1)/2 (om Cosh[2 ombig]+1)]; mu4 = (om-1)^2/8 (om^2 (om^4+2 om^3+3 om^2-3) Cosh[4 ombig]+4 om^2 (om+2) Cosh[2 ombig]+3 (2 om+1)); mu4/sig^4 //N] /; ParameterQ[UnboundedJohnsonDistribution[g,d]]



(* Functionals of Bounded distribution. *)

BoundedJohnsonDistribution/: ParameterQ[BoundedJohnsonDistribution[g_, d_]] := And[ If[FreeQ[N[p], Complex], True, Message[BoundedJohnsoDistribution::realparm, p]; False], If[FreeQ[N[g], Complex], True, Message[BoundedJohnsoDistribution::realparm, d]; False], If[N[d] > 0, True, Message[BoundedJohnsonDistribution::posparm, d]; False, True]]

BoundedJohnsonDistribution::realparm = "Parameter `1` is expected to be real."
BoundedJohnsonDistribution::posparm = "Parameter `1` is expected to be positive."

BoundedJohnsonDistribution/: Domain[BoundedJohnsonDistribution[g_, d_]] := Interval[{0,1}]

BoundedJohnsonDistribution/: DomainQ[BoundedJohnsonDistribution[g_, d_], x_] := IntervalMemberQ[Domain[BoundedJohnsonDistribution[g, d]], x]

BoundedJohnsonDistribution/: PDF[BoundedJohnsonDistribution[g_, d_], x_] := 1./Sqrt[2 Pi] Exp[-(g+d Log[x/(1-x)])^2/2] d/x/(1-x) /; DomainQ[BoundedJohnsonDistribution[g, d], x] && ParameterQ[BoundedJohnsonDistribution[g, d]]

BoundedJohnsonDistribution/: CDF[BoundedJohnsonDistribution[g_, d_], x_] := 0.5 (Erf[(g+d Log[x/(1-x)])/Sqrt[2.]] + 1.0) /; DomainQ[BoundedJohnsonDistribution[g, d], x] && ParameterQ[BoundedJohnsonDistribution[g,d]]

BoundedJohnsonDistribution/: Quantile[BoundedJohnsonDistribution[g_, d_], fraction_] := Module[{z}, z = Exp[(Quantile[NormalDistribution[0, 1], fraction] - g)/d]; z/(1 + z)] /; QuantileQ[fraction] && ParameterQ[BoundedJohnsonDistribution[g,d]]

BoundedJohnsonDistribution/: Random[BoundedJohnsonDistribution[g_, d_]] := Quantile[BoundedJohnsonDistribution[g, d], Random[]]

BoundedJohnsonDistribution/: RandomArray[BoundedJohnsonDistribution[g_, d_], dim_] := Table[Quantile[BoundedJohnsonDistribution[g, d], Random[]], {dim}] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]

BoundedJohnsonDistribution/: Mean[BoundedJohnsonDistribution[g_, d_], opts___?OptionQ] := Module[{x, ret}, With[{optn = FilterOptions[NIntegrate, Sequence@@Flatten[{opts}]]}, ret = NIntegrate[ x PDF[BoundedJohnsonDistribution[g, d], x], {x, 0, 1}, optn]]; Return[ret]] /; ParameterQ[BoundedJohnsonDistribution[g,d]]

(* Mean expressed as infinite series (see Johnson's paper)
Numerical summation here doesn't seem to be efficient...

BoundedJohnsonDistribution/: Mean[BoundedJohnsonDistribution[g_, d_], opts___?OptionQ] := Module[{n, sum1, sum2, sum3}, sum1 = NSum[Exp[-n^2/2/d^2] Cosh[n (1-2 g d)/2/d^2] Sech[n/2/d^2], {n, 1, Infinity}, opts]; sum2 = NSum[Exp[-(2 n-1)^2 Pi^2 d^2/2] Sin[(2 n-1) Pi g d] Csch[(2 n-1) Pi^2 d^2], {n, 1, Infinity}, opts]; sum3 = NSum[Exp[-2 n^2 Pi^2 d^2] Cos[2 n Pi g d], {n, 1, Infinity}, opts]; 1/Sqrt[2 Pi] Exp[-g^2/2] (1/2/d+sum1/d-2 Pi d sum2)/(1 +2 sum3)] /; ParameterQ[BoundedJohnsonDistribution[g,d]]
*)

(* It should be possible to use recurrence relations to calculate higher moments, as instructed in Johnson's paper. *)

BoundedJohnsonDistribution/: Variance[BoundedJohnsonDistribution[g_, d_], opts___?OptionQ] := Module[{mu,x,ret}, mu = Mean[BoundedJohnsonDistribution[g, d], opts]; With[{optn = FilterOptions[NIntegrate, Sequence@@Flatten[{opts}]]}, ret = NIntegrate[ (x-mu)^2 PDF[BoundedJohnsonDistribution[g, d], x], {x, 0, 1}, optn]]; Return[ret]]  /; ParameterQ[BoundedJohnsonDistribution[g,d]]

BoundedJohnsonDistribution/: StandardDeviation[BoundedJohnsonDistribution[g_, d_], opts___?OptionQ] := Module[{mu,x,ret}, mu = Mean[BoundedJohnsonDistribution[g, d], opts]; With[{optn = FilterOptions[NIntegrate, Sequence@@Flatten[{opts}]]}, ret = NIntegrate[ (x-mu)^2 PDF[BoundedJohnsonDistribution[g, d], x], {x, 0, 1}, optn]]; Sqrt[ret]] /; ParameterQ[BoundedJohnsonDistribution[g,d]]

BoundedJohnsonDistribution/: Skewness[BoundedJohnsonDistribution[g_, d_], opts___?OptionQ] := Module[{sig,mu,x,ret}, mu = Mean[BoundedJohnsonDistribution[g, d], opts]; With[{optn = FilterOptions[NIntegrate, Sequence@@Flatten[{opts}]]}, ret = NIntegrate[ (x-mu)^2 PDF[BoundedJohnsonDistribution[g, d], x], {x, 0, 1}, optn]]; sig = Sqrt[ret]; With[{optn = FilterOptions[NIntegrate, Sequence@@Flatten[{opts}]]}, ret = NIntegrate[ (x-mu)^3 PDF[BoundedJohnsonDistribution[g, d], x], {x, 0, 1}, optn]]; 1/sig^3 ret] /; ParameterQ[BoundedJohnsonDistribution[g,d]]

BoundedJohnsonDistribution/: Kurtosis[BoundedJohnsonDistribution[g_, d_], opts___?OptionQ] := Module[{sig,mu,ret}, mu = Mean[BoundedJohnsonDistribution[g, d], opts]; With[{optn = FilterOptions[NIntegrate, Sequence@@Flatten[{opts}]]}, ret = NIntegrate[ (x-mu)^2 PDF[BoundedJohnsonDistribution[g, d], x], {x, 0, 1}, optn]]; sig = Sqrt[ret]; With[{optn = FilterOptions[NIntegrate, Sequence@@Flatten[{opts}]]}, ret = NIntegrate[ (x-mu)^4 PDF[BoundedJohnsonDistribution[g, d], x], {x, 0, 1}, optn]]; 1/sig^4 ret] /; ParameterQ[BoundedJohnsonDistribution[g,d]]



End[]

Protect[UnboundedJohnsonDistribution, BoundedJohnsonDistribution];

EndPackage[]