70  Probability Density Functions (NM742)

Introduction to NONMEM 7.6.0

A series of built in probability density functions are now available. For a given probability density there is also a cumulative distribution function (densityCDF), and random number generating function (density_rng). The following provides an example of how the exponentially modified normal density (EXPMODNORMAL, which returns –log(density)) may be used (other examples are in ..\examples\densities):

ran_expmodnormal.ctl:

$PROB TESTING RANDOM SAMPLERS
$ABBR FUNCTION EXPMODNORMAL(VQI,10)
$ABBR VECTOR VV(10)
$INPUT AMT TVAL DV
$DATA rsampler.csv

$PRED
QM=theta(1)
SIGV=THETA(2)
LAMBDA=THETA(3)
IF(ICALL==4) THEN
; Simulation block. Enter arguments into vector VV().
The first argument usually the random sample to be generated
VV(1)=1.0 ; set to 1 as place holder
VV(2)=QM
VV(3)=SIGV
VV(4)=LAMBDA
; Use normal source 3, and uniform source 2, defined with $SIML
" CALL EXPMODNORMAL_RNG(3,2,VV)
; The generated sample is stored in VV(1), so set DV to this value
Y=VV(1)
DV=Y
ELSE
; Estimation block. The density function returns –log(pdf), so multiply by 2
to make it a -2LL.
VQI(1)=DV ; data value
VQI(2)=QM ; next are the parameters to the density.
VQI(3)=SIGV
VQI(4)=LAMBDA
; density function is called, using VQI vector as first argument.
; NMTRAN adds arguments to capture first and second derivatives.
WW=EXPMODNORMAL(VQI)
Y=2.0*WW
ENDIF

$THETA 30.0 (0.0,5.0) (0.0,0.4)

$SIMULATION (567811 NORMAL) (2933012 UNIFORM) (445678 NORMAL) SUBPROBLEMS=1
$EST METHOD=0 MAXEVAL=9999 PRINT=1 -2LL NOTHETABOUNDTEST
$COVR
$TABLE TVAL DV NOAPPEND NOPRINT FILE=ran_expmodnormal.tab

The density and densityCDF functions have arguments that are compatible with the FUNC system, in which function provides derivatives (XD), and second derivatives (XDD) (see I.78 Expanded Syntax and Capacity for User-Defined Functions (FUNCA) (NM74)). Thus, even random (eta associated) variables may serve as arguments to the parameters of the density functions. The source code of these densities are in ..\source\DISTRIB.f90, DISTRIBCDF.f90, and DISTRIBRNG.f90. Note that multi-variate densities do not have a corresponding CDF routine.

In the above example, the sources for the random number generator are defined with the $SIML record, and use the same generators as the RANDOM() subroutine. The RANDMTU random number generating system (see RANDMTU Function (NM75) in section I.8 Expansions on Abbreviated and Verbatim Code and Other Items (NM72,NM73,NM74,NM75)) may also be used with these probability density function sample generators. The RANDMTU() routine should be first directly called to initialize the source and seed, and subsequently the particular _RNG() subroutine may be called, making the source argument negative (or 0). Using the example of EXPMODNORMAL:

IF(NEWIND==0) THEN
RVAL=RANDMTU(1,2,5678,-1)
RVAL=RANDMTU(0,3,5678,-1)
ENDIF

   VV(1)=1.0
   VV(2)=0.0
   VV(3)=2.0
   VV(4)=3.5
" CALL EXPMODNORMAL_RNG(-2,-3,VV)
   RVAL=VV(1)

The EXPMODNORMAL_RNG function will call source 2 of RANDMTU(), and obtain the appropriate normal sample and uniform sample.

Below are the list of densities (these routines actually return –log(density)), which are modeled after the format from the Stan manual (Stan Development Team 2015).

1 BERNOULLI

Given \(\theta \in [0,1]\), then for \(y=\{0,1\}\),

\[\text{Bernoulli}(y \mid \theta) = \begin{cases} \theta & \text{if } y=1 \\ 1-\theta & \text{if } y=0 \end{cases}\]

2 BERNOULLILOGIT

Given \(\alpha \in R\), then for \(y=\{0,1\}\),

\[ \text{BernoulliLogit}(y \mid \alpha) = \begin{cases} \log\text{it}^{-1}(\alpha) & \text{if } y=1 \\ 1 - \log\text{it}^{-1}(\alpha) & \text{if } y=0 \end{cases} \]

\[\text{logit}^{-1}(\alpha) = \frac{\exp(\alpha)}{1+\exp(\alpha)}\]

3 BINOMIAL

Given integer \(N\), \(\theta \in [0,1]\), and \(n=\{0,\ldots,N\}\)

\[\text{Binomial}(n \mid N, \theta) = \binom{N}{n} \theta^n (1-\theta)^{N-n}\]

4 BINOMIALLOGIT

Given integer \(N\), \(\alpha \in R\), and \(n=\{0,\ldots,N\}\)

\[\text{Binomial logit}(n \mid N, \alpha) = \binom{N}{n} \text{logit}^{-1}(\alpha)^n (1-\text{logit}^{-1}(\alpha))^{N-n}\]

5 BETABINOMIAL

Given integer \(N\), \(\alpha \in R^+\), \(\beta \in R^+\), and \(n=\{0,\ldots,N\}\)

\[\text{Beta Binomial}(n \mid N, \alpha, \beta) = \binom{N}{n} \frac{B(n+\alpha, N-n+\beta)}{B(\alpha, \beta)}\]

\[B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}\]

\[\Gamma(x) = (x+1)!\]

6 HYPERGEOMETRIC

Given integer \(N\), integer \(a\), integer \(b\), and \(n=\{\max(0,N-b),\ldots,\min(a,N)\}\),

\[\text{Hypergeometric}(n \mid N, a, b) = \frac{\binom{a}{n}\binom{b}{N-n}}{\binom{a+b}{N}}\]

7 CATEGORICAL

For positive integer \(N\), \(\theta \in R^N\) are \(N\) non-negative set of values summing to one (\(N\)-simplex), then for \(y=\{1,\ldots,N\}\),

\[\text{Categorical}(y \mid \theta) = \theta_y\]

Values are entered as:

X(1)=y
X(2)=N
X(3)=th(1)
X(4)=th(2)
…
X(N+1)=th(N-1)

8 CATEGORICALLOGIT

For positive integer \(N\), \(\beta \in R^N\), then for \(y=\{1,\ldots,N\}\),

\[\theta_y = \frac{\exp(\beta_y)}{\sum_{k=1}^{N} \exp(\beta_k)}\]

\[\text{Categorical}(y \mid \theta) = \theta_y\]

Values are entered as:

X(1)=y
X(2)=N
X(3)=beta(1)
X(4)=beta(2)
…
X(N+1)=beta(N-1)

and Beta(N) is assumed 0.

9 ORDEREDLOGISTIC

For integer \(K\), such that \(c_k<c_{k+1}\) for \(k=\{1,\ldots,K-2\}\), and \(\eta \in R\), then for \(k=\{1,\ldots,K\}\)

\[ \text{OrderedLogistic}(k \mid \eta, c) = \begin{cases} 1 - \text{logit}^{-1}(\eta - c_1) & \text{if } k=1 \\ \text{logit}^{-1}(\eta - c_{k-1}) - \text{logit}^{-1}(\eta - c_k) & \text{if } 1 < k < K \\ \text{logit}^{-1}(\eta - c_{k-1}) & \text{if } k=K \end{cases} \]

10 NEGBINOMIAL

For \(\alpha \in R^+\), and \(\beta \in R^+\), then for positive integer \(y\)

\[ \text{NegBinomial}(y \mid \alpha, \beta) = \binom{y+\alpha-1}{\alpha-1} \left(\frac{\beta}{\beta+1}\right)^{\alpha} \left(\frac{1}{\beta+1}\right)^{y} \]

11 NEGBINOMIAL2

For \(\mu \in R^+\), and \(\phi \in R^+\), then for positive integer \(y\)

\[ \text{NegBinomial2}(y \mid \mu, \phi) = \binom{y+\phi-1}{\phi-1} \left(\frac{\mu}{\mu+\phi}\right)^{y} \left(\frac{\phi}{\mu+\phi}\right)^{\phi} \]

12 NEGBINOMIAL2LOG

For \(\eta = \log(\mu) \in R\), and \(\phi \in R^+\), then for integer positive \(y\)

\[\text{NegBinomial2log}(y \mid \eta, \phi) = \text{NegBinomial2}(y \mid \exp(\exp(\eta)), \phi)\]

13 POISSON

For \(\lambda \in R^+\), then for integer \(n \geq 0\)

\[\text{Poisson}(n \mid \lambda) = \frac{1}{n!} \lambda^n \exp(-\lambda)\]

14 POISSONLOG

For \(\alpha = \log(\lambda) \in R\), then for integer \(n \geq 0\)

\[\text{Poisson}(n \mid \alpha) = \frac{1}{n!} \exp(n\alpha - \exp(\alpha))\]

15 MULTINOMIAL

If \(K\) is positive integer, \(\theta \in K-\text{simplex}\), then for integer \(y\), where

\[\sum_{k=1}^{K} y_k = N\]

\[\text{Multinomial}(y \mid \theta) = \binom{N}{y_1,\ldots,y_k} \prod_{k=1}^{K} \theta_k^{y_k}\]

Enter values in X array as follows for the –log density multinomial:

X(1)=K
X(2)=y(1)
…
X(K+1)=y(K)
X(K+2)=TH(1)
..
X(2K)=TH(K-1)

For the random number generator,

X(1)=K
X(2)=N ; Enter total value
X(3)=1.0 ; place holder (multinomial_rng will fill in X(3) to X(K+1)
…
X(K+1)=1.0
X(K+2)=TH(1)
..
X(2K)=TH(K-1)

As a multivariate density, Multinomial does not have a CDF routine.

16 NORMAL

For \(\mu \in R\), \(\sigma \in R^+\), \(y \in R\),

\[ \text{NORMAL}(y \mid \mu, \sigma) = \frac{1}{\sqrt{2\pi}\sigma} \exp\left(-\frac{1}{2}\left(\frac{y-\mu}{\sigma}\right)^2\right) \]

17 EXPMODNORMAL

For \(\mu \in R\), \(\sigma \in R^+\), \(\lambda \in R^+\), \(y \in R\),

\[ \text{EXPMODNORMAL}(y \mid \mu, \sigma, \lambda) = \frac{\lambda}{2} \exp\left(\frac{\lambda}{2}(2\mu + \lambda\sigma^2 - 2y)\right) \cdot 2\left(1 - \Phi\left(\frac{\mu + \lambda\sigma^2 - y}{\sigma}\right)\right) \]

18 SKEWNORMAL

For \(\mu \in R\), \(\sigma \in R^+\), \(\alpha \in R\), \(y \in R\),

\[ \text{SKEWNORMAL}(y \mid \mu, \sigma, \alpha) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left(-\frac{1}{2}\left(\frac{y-\mu}{\sigma}\right)^2\right) \cdot 2\Phi\left(\alpha \frac{y-\mu}{\sigma}\right) \]

19 STUDENTT

For \(\mu \in R\), \(\nu \in R^+\), \(\sigma \in R^+\), \(y \in R\),

\[ \text{STUDENTT}(y \mid \nu, \mu, \sigma) = \frac{\Gamma((\nu+1)/2)}{\Gamma(\nu/2)} \frac{1}{\sqrt{\nu\pi}\sigma} \left(1 + \frac{1}{\nu}\left(\frac{y-\mu}{\sigma}\right)^2\right)^{-(\nu+1)/2} \]

Likelihood estimator STUDENTT2 in DISTRIB.f90 is an alternative coding, but otherwise equivalent to STUDENTT.

There are three types of STUDENTT random number generators available in DISTRIBRNG.f90:

STUDDENTT_RNG(K,X): Given uniform Random generator K, and parameters X(2)..., return random variable X(1). Uses the TDEV2 routine in GENERAL.f90, which uses efficient random number generators for NU=1, 2, 4, 6, and 10. NU=X(2) must be integer valued.

STUDDENTT2_RNG(K,X): Given Normal Random generator K, and parameters X(2)..., return random variable X(1). NU=X(2) may be non-integer. Algorithm uses two normal random variables to generate one univariate t-variable.

STUDDENTT3_RNG(K,X): Given normal generator uniform Random generator K1, uniform generator K2, and parameters X(2)..., return random variable X(1). NU=X(2) may by non-integer. This routine uses the normal variable/square root of chi-square variable) algorithm.

There are also the STUDENTTB series routines in DISTRIB.f90 and DISTRIB_RNG.f90, which deal with pairs of random correlated t-distributed samples. Please read comments in the source code for more information.

20 DOUBLEEXPONENTIAL

For \(\mu \in R\), \(\sigma \in R^+\), \(y \in R\),

\[\text{DOUBLEEXPONENTIAL}(y \mid \mu, \sigma) = \frac{1}{2\sigma} \exp\left(-\frac{|y-\mu|}{\sigma}\right)\]

21 LOGISTIC

For \(\mu \in R\), \(\sigma \in R^+\), \(y \in R\),

\[ \text{LOGISTIC}(y \mid \mu, \sigma) = \frac{1}{\sigma} \exp\left(-\frac{y-\mu}{\sigma}\right) \left(1 + \exp\left(-\frac{y-\mu}{\sigma}\right)\right)^{-2} \]

22 GUMBEL

For \(\mu \in R\), \(\beta \in R^+\), \(y \in R\),

\[\text{GUMBEL}(y \mid \mu, \beta) = \frac{1}{\beta} \exp\left(-\frac{y-\mu}{\beta} - \exp\left(-\frac{y-\mu}{\beta}\right)\right)\]

23 LOGNORMAL

For \(\mu \in R\), \(\sigma \in R^+\), \(y \in R\),

\[\text{LOGNORMAL}(y \mid \mu, \sigma) = \frac{1}{\sqrt{2\pi}\sigma y} \exp\left(-\frac{1}{2}\left(\frac{\log(y)-\mu}{\sigma}\right)^2\right)\]

24 CHISQUARE

For \(\nu \in R\), \(y \in R^+\),

\[\text{CHISQUARE}(y \mid \nu) = \frac{2^{-\nu/2}}{\Gamma(\nu/2)} y^{\nu/2-1} \exp\left(-\frac{1}{2}y\right)\]

25 INVCHISQUARE

For \(\nu \in R\), \(y \in R^+\),

\[\text{INVCHISQUARE}(y \mid \nu) = \frac{2^{-\nu/2}}{\Gamma(\nu/2)} y^{-\nu/2+1} \exp\left(-\frac{1}{2y}\right)\]

26 SCALEDINVCHISQUARE

For \(\nu \in R\), \(y \in R^+\),

\[ \text{SCALEDINVCHISQUARE}(y \mid \nu, \sigma) = \frac{(\nu/2)^{\nu/2}}{\Gamma(\nu/2)} \sigma^{\nu} y^{-\nu/2-1} \exp\left(-\frac{1}{2}\nu\sigma^2 \frac{1}{y}\right) \]

27 EXPONENTIAL

For \(\beta \in R^+\), \(y \in R\),

\[\text{EXPONENTIAL}(y \mid \beta) = \beta \exp(-\beta y)\]

28 GAMMA

For \(\alpha \in R^+\), \(\beta \in R^+\), \(y \in R^+\),

\[\text{GAMMA}(y \mid \alpha, \beta) = \frac{\beta^{\alpha}}{\Gamma(\alpha)} y^{\alpha-1} \exp(-\beta y)\]

When using gfortran compiler, you will need to add

" first
" external gamma

in the first lines of $PK or $PRED, so the external gamma function supplied in NONMEM is used and not the intrinsic gamma function.

29 INVGAMMA

For \(\alpha \in R^+\), \(\beta \in R^+\), \(y \in R^+\),

\[\text{INVGAMMA}(y \mid \alpha, \beta) = \frac{\beta^{\alpha}}{\Gamma(\alpha)} y^{-\alpha-1} \exp\left(-\beta \frac{1}{y}\right)\]

30 WEIBULL

For \(\alpha \in R^+\), \(\sigma \in R^+\), \(y \in R^+\),

\[ \text{WEIBULL}(y \mid \alpha, \sigma) = \frac{\alpha}{\sigma} \left(\frac{y}{\sigma}\right)^{\alpha-1} \exp\left(-\left(\frac{y}{\sigma}\right)^{\alpha}\right) \]

31 FRECHET

For \(\alpha \in R^+\), \(\sigma \in R^+\), \(y \in R^+\),

\[ \text{FRECHET}(y \mid \alpha, \sigma) = \frac{\alpha}{\sigma} \left(\frac{y}{\sigma}\right)^{-\alpha-1} \exp\left(-\left(\frac{y}{\sigma}\right)^{-\alpha}\right) \]

32 RAYLEIGH

For \(\sigma \in R^+\), \(y \in R^+\),

\[\text{RAYLEIGH}(y \mid \sigma) = \frac{y}{\sigma^2} \exp\left(-\frac{y^2}{2\sigma^2}\right)\]

33 PARETO

For real positive \(y_{\min}\), \(\alpha \in R^+\), \(y \geq y_{\min}\),

\[\text{PARETO}(y \mid y_{\min}, \alpha) = \frac{\alpha y_{\min}^{\alpha}}{y^{\alpha+1}}\]

34 PARETO2

For \(\mu \in R\), \(\lambda \in R^+\), \(\alpha \in R^+\), \(y \geq \mu\),

\[\text{PARETO2}(y \mid \mu, \lambda, \alpha) = \frac{\alpha}{\lambda} \left(1 + \frac{y-\mu}{\lambda}\right)^{-(\alpha+1)}\]

35 BETA

For \(\alpha \in R^+\), \(\beta \in R^+\), \(\theta \in (0,1)\)

\[\text{BETA}(\theta \mid \alpha, \beta) = \frac{1}{B(\alpha, \beta)} \theta^{\alpha-1} (1-\theta)^{\beta-1}\]

36 DIRICHLET

For positive integer \(K\), \(\alpha \in (R^+)^K\), then for \(\theta \in K-\text{simplex}\)

\[ \text{DIRICHLET}(\theta \mid \alpha) = \frac{ \Gamma\left(\sum_{k=1}^{K} \alpha_k\right) }{ \prod_{k=1}^{K} \Gamma(\alpha_k) } \prod_{k=1}^{K} \theta_k^{\alpha_k-1} \]

Enter values as follows:

X(1)=K
X(2)=TH(1)
…
X(K)=TH(K-1)
X(K+1)=alpha(1)
…
X(2K)=alpha(K)

As a multivariate density, DIRICHLET does not have a CDF routine.

37 VON MISES

For positive real \(\mu\), real positive \(\kappa\), and real \(-\pi \leq y-\mu \leq \pi\),

\[\text{VonMises}(y \mid \mu, \kappa) = \frac{\exp(\kappa \cos(y-\mu))}{2\pi I_0(\kappa)}\]

Where \(I_0(\kappa)\) is a modified Bessel function. There is no VonMisesCDF routine available. Two methods of random sample generation are available:

Vonmises_rng: Lucio Barabesi, Statistica Applicata, Vol 7 number 4, 1995, pp. 417-426, page 420,

Vonmises_rng2: Stan system (Stan Development Team 2015), which in turn is based on D. J. Best and N. I. Fisher, Journal of the Royal Statistical Society. Series C (Applied Statistics), Vol. 28, No. 2 (1979), pp. 152-157

Stan Development Team. 2015. Stan Modeling Language Users Guide and Reference Manual, Version 2.6.3. https://mc-stan.org/docs/.