/*********  Laplace-Distributed Random Numbers ***********
Author: Mico Loretan  mailto:loretanm@frb.gov
Date: 3 Aug 1999
Caveat: NO WARRANTIES WHATSOEVER, either explicit or implicit.

This code generates Laplace-Distributed Random Numbers
using the inverse-cdf technique,
I'm enclosing a sample proc that performs this operation in
Gauss.  Feel free to redistribute/critique/improve the code.

***********************************************************/

/* laplace.src */



proc laplace(m,n,a,b);
/*
:: Mico Loretan, August 1999.
::
:: This proc returns an m x n matrix of Laplace-distributed
:: random variables, with parameters "a" and "b".
:: We have: E(X) = a and Var(X)=2*b^2.
::
:: Inputs:  m   number of rows of output matrix
::  n   number of columns of output matrix
::  (m and n must be positive integers)
::  a   Location parameter of Laplace r.v.
::  (ExE-conformable with m and n)
::  b   Dispersion parameter of Laplace r.v.
::  (ExE-conformable with m and n)
:: Output:  la  mxn matrix of Laplace-distributed r.v.'s
::
:: Note:If r.v.'s with zero mean and standard scale
::  (i.e., Var(x)=1) are desired, set a=0 and
::  b=sqrt(1/2).
*/

local u,u1,u2, e, la;

if (b<=0);
errorlog "Invalid input in proc laplace() -- nonpositive
dispersion!";
end;
endif;

u = rndu(m,n);  /* mxn matrix of standard uniform r.v.'s */
e = u .< 0.5;   /* mxn matrix of indicator r.v.'s:
   e_[ij]=1 if u_[ij]<0.5 */

u1 = missex(u,1-e);   /* u1_[ij] = missing value if u_[ij]>0.5 */
u2 = missex(u,e); /* u2_[ij] = missing value if u_[ij]<0.5 */

disable;/* suppress complaints about ops with missing values */
la = missrv(ln(2*u1),0) - missrv(ln(2*(1-u2)),0);
enable;
ndpclex;

la = a+b.*la;  /* finally, adjust for location and dispersion */
retp(la);
endp; /* laplace() */