1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: c_laplace.c *
8 * *
9 * REFERENCES: *
10 * *
11 * [3] N.L. Johnson, S. Kotz and N. Balakrishnan *
12 * Continuous Univariate Distributions, *
13 * Volume 2, 2nd edition *
14 * John Wiley & Sons, Inc., New York, 1995 *
15 * *
16 *****************************************************************************
17 * *
18 * Laplace distribution [3; ch.24, p.164] *
19 * *
20 * pdf: f(x) = exp(- abs(x-theta) / phi ) *
21 * domain: -infinity < x < infinity *
22 * constant: 1 / (2 * phi) *
23 * *
24 * parameters: *
25 * 0: theta ... location *
26 * 1: phi > 0 ... scale *
27 * *
28 *****************************************************************************
29 * *
30 * standard form *
31 * *
32 * pdf: f(x) = exp(- abs(x)) *
33 * domain: -infinity < x < infinity *
34 * constant: 2 *
35 * *
36 * parameters: none *
37 * *
38 * 0: theta = 0 ... location *
39 * 1: phi = 1 ... scale *
40 * *
41 *****************************************************************************
42 * *
43 * Copyright (c) 2000-2006 Wolfgang Hoermann and Josef Leydold *
44 * Department of Statistics and Mathematics, WU Wien, Austria *
45 * *
46 * This program is free software; you can redistribute it and/or modify *
47 * it under the terms of the GNU General Public License as published by *
48 * the Free Software Foundation; either version 2 of the License, or *
49 * (at your option) any later version. *
50 * *
51 * This program is distributed in the hope that it will be useful, *
52 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
53 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
54 * GNU General Public License for more details. *
55 * *
56 * You should have received a copy of the GNU General Public License *
57 * along with this program; if not, write to the *
58 * Free Software Foundation, Inc., *
59 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
60 * *
61 *****************************************************************************/
62
63 /*---------------------------------------------------------------------------*/
64
65 #include <unur_source.h>
66 #include <distr/distr_source.h>
67 #include <distr/cont.h>
68 #include <specfunct/unur_specfunct_source.h>
69 #include "unur_distributions.h"
70 #include "unur_distributions_source.h"
71 #include "unur_stddistr.h"
72
73 /*---------------------------------------------------------------------------*/
74 static const char distr_name[] = "laplace";
75
76 /* parameters */
77 #define theta params[0]
78 #define phi params[1]
79
80 #define DISTR distr->data.cont
81 /* #define NORMCONSTANT (distr->data.cont.norm_constant) */
82
83 /* function prototypes */
84 static double _unur_pdf_laplace( double x, const UNUR_DISTR *distr );
85 static double _unur_dpdf_laplace( double x, const UNUR_DISTR *distr );
86 static double _unur_cdf_laplace( double x, const UNUR_DISTR *distr );
87 static double _unur_invcdf_laplace( double u, const UNUR_DISTR *distr );
88 static double _unur_logpdf_laplace( double x, const UNUR_DISTR *distr );
89 static double _unur_dlogpdf_laplace( double x, const UNUR_DISTR *distr );
90
91 static int _unur_upd_mode_laplace( UNUR_DISTR *distr );
92 static int _unur_upd_area_laplace( UNUR_DISTR *distr );
93 static int _unur_set_params_laplace( UNUR_DISTR *distr, const double *params, int n_params );
94
95 /*---------------------------------------------------------------------------*/
96
97 double
_unur_pdf_laplace(double x,const UNUR_DISTR * distr)98 _unur_pdf_laplace( double x, const UNUR_DISTR *distr )
99 {
100 register const double *params = DISTR.params;
101
102 if (x<theta)
103 return ( exp( (x-theta)/phi ) / (2.*phi) );
104 /* else */
105 return ( exp( (theta-x)/phi ) / (2.*phi) );
106
107 } /* end of _unur_pdf_laplace() */
108
109 /*---------------------------------------------------------------------------*/
110
111 double
_unur_logpdf_laplace(double x,const UNUR_DISTR * distr)112 _unur_logpdf_laplace( double x, const UNUR_DISTR *distr )
113 {
114 register const double *params = DISTR.params;
115
116 if (x<theta)
117 return ( (x-theta)/phi - log(2.*phi) );
118 /* else */
119 return ( (theta-x)/phi - log(2.*phi) );
120
121 } /* end of _unur_logpdf_laplace() */
122
123 /*---------------------------------------------------------------------------*/
124
125 double
_unur_dpdf_laplace(double x,const UNUR_DISTR * distr)126 _unur_dpdf_laplace( double x, const UNUR_DISTR *distr )
127 {
128 register const double *params = DISTR.params;
129 register double z;
130
131 z = (x>theta) ? (x-theta)/phi : (theta-x)/phi;
132
133 if (_unur_iszero(z)) /* derivative is not defined, but ... */
134 return 0.; /* a tangent parallel to x-axis is possible. */
135
136 return ( ((x>theta) ? -exp(-z)/phi : exp(-z)/phi) / (2.*phi) );
137 } /* end of unur_dpdf_laplace() */
138
139 /*---------------------------------------------------------------------------*/
140
141 double
_unur_dlogpdf_laplace(double x,const UNUR_DISTR * distr)142 _unur_dlogpdf_laplace( double x, const UNUR_DISTR *distr )
143 {
144 register const double *params = DISTR.params;
145 return ( (x<theta) ? 1./phi : -1./phi );
146 } /* end of unur_dlogpdf_laplace() */
147
148 /*---------------------------------------------------------------------------*/
149
150 double
_unur_cdf_laplace(double x,const UNUR_DISTR * distr)151 _unur_cdf_laplace( double x, const UNUR_DISTR *distr )
152 {
153 register const double *params = DISTR.params;
154 register double z;
155
156 z = (x-theta)/phi;
157 return ( (x>theta) ? 1.-0.5 * exp(-z) : 0.5*exp(z) );
158 } /* end of _unur_cdf_laplace() */
159
160 /*---------------------------------------------------------------------------*/
161
162 double
_unur_invcdf_laplace(double U,const UNUR_DISTR * distr)163 _unur_invcdf_laplace( double U, const UNUR_DISTR *distr )
164 {
165 register const double *params = DISTR.params;
166 double X;
167
168 U *= 2.;
169 X = (U>1.) ? -log(2.-U) : log(U);
170 return ((DISTR.n_params==0) ? X : theta + phi * X );
171 } /* end of _unur_invcdf_laplace() */
172
173 /*---------------------------------------------------------------------------*/
174
175 int
_unur_upd_mode_laplace(UNUR_DISTR * distr)176 _unur_upd_mode_laplace( UNUR_DISTR *distr )
177 {
178 DISTR.mode = DISTR.theta;
179
180 /* mode must be in domain */
181 if (DISTR.mode < DISTR.domain[0])
182 DISTR.mode = DISTR.domain[0];
183 else if (DISTR.mode > DISTR.domain[1])
184 DISTR.mode = DISTR.domain[1];
185
186 return UNUR_SUCCESS;
187 } /* end of _unur_upd_mode_laplace() */
188
189 /*---------------------------------------------------------------------------*/
190
191 int
_unur_upd_area_laplace(UNUR_DISTR * distr)192 _unur_upd_area_laplace( UNUR_DISTR *distr )
193 {
194 /* normalization constant: none */
195
196 if (distr->set & UNUR_DISTR_SET_STDDOMAIN) {
197 DISTR.area = 1.;
198 return UNUR_SUCCESS;
199 }
200
201 /* else */
202 DISTR.area = ( _unur_cdf_laplace( DISTR.domain[1],distr)
203 - _unur_cdf_laplace( DISTR.domain[0],distr) );
204 return UNUR_SUCCESS;
205
206 } /* end of _unur_upd_area_laplace() */
207
208 /*---------------------------------------------------------------------------*/
209
210 int
_unur_set_params_laplace(UNUR_DISTR * distr,const double * params,int n_params)211 _unur_set_params_laplace( UNUR_DISTR *distr, const double *params, int n_params )
212 {
213 /* check number of parameters for distribution */
214 if (n_params < 0) n_params = 0;
215 if (n_params > 2) {
216 _unur_warning(distr_name,UNUR_ERR_DISTR_NPARAMS,"too many");
217 n_params = 2; }
218 if (n_params > 0)
219 CHECK_NULL(params,UNUR_ERR_NULL);
220
221 /* check parameter phi */
222 if (n_params == 2 && phi <= 0.) {
223 _unur_error(distr_name,UNUR_ERR_DISTR_DOMAIN,"phi <= 0");
224 return UNUR_ERR_DISTR_DOMAIN;
225 }
226
227 /* copy parameters for standard form: none */
228
229 /* default parameters */
230 DISTR.theta = 0.;
231 DISTR.phi = 1.;
232
233 /* copy optional parameters */
234 switch (n_params) {
235 case 2:
236 DISTR.phi = phi;
237 case 1:
238 DISTR.theta = theta;
239 default:
240 n_params = 2; /* number of parameters for non-standard form */
241 }
242
243 /* store number of parameters */
244 DISTR.n_params = n_params;
245
246 /* set (standard) domain */
247 if (distr->set & UNUR_DISTR_SET_STDDOMAIN) {
248 DISTR.domain[0] = -INFINITY; /* left boundary */
249 DISTR.domain[1] = INFINITY; /* right boundary */
250 }
251
252 return UNUR_SUCCESS;
253 } /* end of _unur_set_params_laplace() */
254
255 /*---------------------------------------------------------------------------*/
256
257 struct unur_distr *
unur_distr_laplace(const double * params,int n_params)258 unur_distr_laplace( const double *params, int n_params )
259 {
260 register struct unur_distr *distr;
261
262 /* get new (empty) distribution object */
263 distr = unur_distr_cont_new();
264
265 /* set distribution id */
266 distr->id = UNUR_DISTR_LAPLACE;
267
268 /* name of distribution */
269 distr->name = distr_name;
270
271 /* how to get special generators */
272 /* DISTR.init = _unur_stdgen_laplace_init; */
273
274 /* functions */
275 DISTR.pdf = _unur_pdf_laplace; /* pointer to PDF */
276 DISTR.logpdf = _unur_logpdf_laplace; /* pointer to logPDF */
277 DISTR.dpdf = _unur_dpdf_laplace; /* pointer to derivative of PDF */
278 DISTR.dlogpdf = _unur_dlogpdf_laplace; /* pointer to deriv. of logPDF */
279 DISTR.cdf = _unur_cdf_laplace; /* pointer to CDF */
280 DISTR.invcdf = _unur_invcdf_laplace; /* pointer to inverse CDF */
281
282 /* indicate which parameters are set */
283 distr->set = ( UNUR_DISTR_SET_DOMAIN |
284 UNUR_DISTR_SET_STDDOMAIN |
285 UNUR_DISTR_SET_PDFAREA |
286 UNUR_DISTR_SET_MODE );
287
288 /* set parameters for distribution */
289 if (_unur_set_params_laplace(distr,params,n_params)!=UNUR_SUCCESS) {
290 free(distr);
291 return NULL;
292 }
293
294 /* normalization constant: none */
295
296 /* mode and area below p.d.f. */
297 DISTR.mode = DISTR.theta;
298 DISTR.area = 1.;
299
300 /* function for setting parameters and updating domain */
301 DISTR.set_params = _unur_set_params_laplace;
302
303 /* function for updating derived parameters */
304 DISTR.upd_mode = _unur_upd_mode_laplace; /* funct for computing mode */
305 DISTR.upd_area = _unur_upd_area_laplace; /* funct for computing area */
306
307 /* return pointer to object */
308 return distr;
309
310 } /* end of unur_distr_laplace() */
311
312 /*---------------------------------------------------------------------------*/
313 #undef theta
314 #undef phi
315 #undef DISTR
316 /*---------------------------------------------------------------------------*/
317