*! Date : 03 Jan 2019 *! Version : 1.4 *! Authors : Michael J Grayling & Adrian P Mander /* 16/04/15 v1.0 Basic version complete. 02/10/15 v1.1 Changed method for computing density for better stability. 08/09/17 v1.2 Added logdensity option and support for non-integer df. 09/09/17 v1.3 Changed so df = 0 or df = . corresponds to mvnormalden. 03/01/19 v1.4 Converted delta and sigma to be optional with internal defaults. */ program define mvtden, rclass version 15.0 syntax , x(numlist) [DELta(numlist) Sigma(string) df(real 1) LOGdensity] ///// Check input variables //////////////////////////////////////////////////// local lenx:list sizeof x local lendelta:list sizeof delta if ((`lendelta' != 0) & (`lenx' != `lendelta')) { di "{error}Vector of quantiles (x) and vector of non-centrality parameters (delta) must be of equal length." exit(198) } if ("`sigma'" != "") { if (colsof(`sigma') != rowsof(`sigma')) { di "{error}Scale matrix Sigma (sigma) must be square." exit(198) } if (`lenx' != colsof(`sigma')) { di "{error}Vector of quantiles (x) must be the same length as the dimension of scale matrix Sigma (sigma)." exit(198) } cap mat chol = cholesky(`sigma') if (_rc > 0) { di "{error}Covariance matrix Sigma (sigma) must be symmetric positive-definite." exit(198) } } else { mat chol = I(`lenx') } if (`df' < 0) { di "{error}Degrees of freedom (df) must be greater than or equal to zero." exit(198) } ///// Perform main computations //////////////////////////////////////////////// local matax "" foreach l of local x { if "`matax'" == "" local matax "`l'" else local matax "`matax',`l'" } mat x = (`matax') if (`lendelta' != 0) { local matadelta "" foreach l of local delta { if "`matadelta'" == "" local matadelta "`l'" else local matadelta "`matadelta',`l'" } mat delta = (`matadelta') } else { mat delta = J(1, `lenx', 0) } mata: mvtden_void(`df') ///// Output /////////////////////////////////////////////////////////////////// return scalar log_density = ret_density[2, 1] return scalar density = ret_density[1, 1] if ("`logdensity'" == "") { di "{txt}density = {res}" ret_density[1, 1] "{txt}" } else { di "{txt}log(density) = {res}" ret_density[2, 1] "{txt}" } end ///// Mata ///////////////////////////////////////////////////////////////////// mata: void mvtden_void(df) { x = st_matrix("x") delta = st_matrix("delta") C = st_matrix("chol") st_matrix("ret_density", mvtden_mata(x, delta, ., df, C)) } real colvector mvtden_mata(real vector x, real vector delta, real matrix Sigma, real scalar df,| real matrix C) { if (args() < 5) { C = cholesky(Sigma) } if ((df == 0) | (df == .)) { log_density = -sum(log(diagonal(C))) - 0.5*(rows(C)*log(2*pi()) + colsum(lusolve(C, (x - delta)'):^2)) } else { k = rows(C) log_density = lngamma((df + k)/2) - (lngamma(df/2) + sum(log(diagonal(C))) + (k/2)*log(pi()*df)) - 0.5*(df + k)* log(1 + (colsum(lusolve(C, (x - delta)'):^2)/df)) } return((exp(log_density) \ log_density)) } end