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