PageRenderTime 57ms CodeModel.GetById 22ms RepoModel.GetById 1ms app.codeStats 0ms

/unuran-1.8.0/src/distributions/c_F.c

#
C | 372 lines | 202 code | 75 blank | 95 comment | 39 complexity | b21d7ae48517a0f8ac82f88cd32fe590 MD5 | raw file
Possible License(s): GPL-2.0
  1. /*****************************************************************************
  2. * *
  3. * UNURAN -- Universal Non-Uniform Random number generator *
  4. * *
  5. *****************************************************************************
  6. * *
  7. * FILE: c_F.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. * distr: F distribution [3; ch.27, p.332] *
  19. * *
  20. * pdf: f(x) = (x^{{nu_1}/2-1}) / ((1+nu_1/nu_2 x)^{(nu_1+nu_2)/2}) *
  21. * domain: 0 < x < infinity *
  22. * constant: \frac{(nu_1/nu_2)^{{nu_1}/2}/Beta({nu_1}/2,{nu_2}/2) *
  23. * *
  24. * parameters: *
  25. * 0: nu_1 > 0 ... shape (degrees of freedom) *
  26. * 1: nu_2 > 0 ... shape (degrees of freedom) *
  27. * *
  28. *****************************************************************************
  29. * *
  30. * Copyright (c) 2000-2010 Wolfgang Hoermann and Josef Leydold *
  31. * Department of Statistics and Mathematics, WU Wien, Austria *
  32. * *
  33. * This program is free software; you can redistribute it and/or modify *
  34. * it under the terms of the GNU General Public License as published by *
  35. * the Free Software Foundation; either version 2 of the License, or *
  36. * (at your option) any later version. *
  37. * *
  38. * This program is distributed in the hope that it will be useful, *
  39. * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  40. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  41. * GNU General Public License for more details. *
  42. * *
  43. * You should have received a copy of the GNU General Public License *
  44. * along with this program; if not, write to the *
  45. * Free Software Foundation, Inc., *
  46. * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
  47. * *
  48. *****************************************************************************/
  49. /*---------------------------------------------------------------------------*/
  50. #include <unur_source.h>
  51. #include <distr/distr_source.h>
  52. #include <distr/cont.h>
  53. #include <specfunct/unur_specfunct_source.h>
  54. #include "unur_distributions.h"
  55. #include "unur_distributions_source.h"
  56. #include "unur_stddistr.h"
  57. /*---------------------------------------------------------------------------*/
  58. static const char distr_name[] = "F";
  59. /*---------------------------------------------------------------------------*/
  60. /* parameters */
  61. #define nua params[0]
  62. #define nub params[1]
  63. /*---------------------------------------------------------------------------*/
  64. #define DISTR distr->data.cont
  65. #define LOGNORMCONSTANT (distr->data.cont.norm_constant)
  66. /*---------------------------------------------------------------------------*/
  67. /* function prototypes */
  68. static double _unur_pdf_F( double x, const UNUR_DISTR *distr );
  69. static double _unur_logpdf_F( double x, const UNUR_DISTR *distr );
  70. static double _unur_dpdf_F( double x, const UNUR_DISTR *distr );
  71. static double _unur_dlogpdf_F( double x, const UNUR_DISTR *distr );
  72. static double _unur_cdf_F( double x, const UNUR_DISTR *distr );
  73. #ifdef _unur_SF_invcdf_F
  74. static double _unur_invcdf_F( double x, const UNUR_DISTR *distr );
  75. #endif
  76. static int _unur_upd_mode_F( UNUR_DISTR *distr );
  77. static int _unur_upd_area_F( UNUR_DISTR *distr );
  78. inline static double _unur_lognormconstant_F( const double *params, int n_params );
  79. static int _unur_set_params_F( UNUR_DISTR *distr, const double *params, int n_params );
  80. /*---------------------------------------------------------------------------*/
  81. double
  82. _unur_pdf_F(double x, const UNUR_DISTR *distr)
  83. {
  84. register const double *params = DISTR.params;
  85. if (x < 0.)
  86. /* out of support */
  87. return 0.;
  88. else if (_unur_iszero(x)) {
  89. if (nua < 2.)
  90. return INFINITY;
  91. else if (_unur_isfsame(nua,2.))
  92. return exp(-LOGNORMCONSTANT);
  93. else
  94. return 0.;
  95. }
  96. else
  97. return exp ((nua/2. - 1.)*log(x) - 0.5*(nua + nub)*log(1. + x * nua / nub) - LOGNORMCONSTANT);
  98. } /* end of _unur_pdf_F() */
  99. /*---------------------------------------------------------------------------*/
  100. double
  101. _unur_logpdf_F(double x, const UNUR_DISTR *distr)
  102. {
  103. register const double *params = DISTR.params;
  104. if (x < 0.)
  105. /* out of support */
  106. return -INFINITY;
  107. else if (_unur_iszero(x)) {
  108. if (nua < 2.)
  109. return INFINITY;
  110. else if (_unur_isfsame(nub,2.))
  111. return -LOGNORMCONSTANT;
  112. else
  113. return -INFINITY;
  114. }
  115. else
  116. return ((nua/2. - 1.)*log(x) - 0.5*(nua + nub)*log(1. + x * nua / nub) - LOGNORMCONSTANT);
  117. } /* end of _unur_logpdf_F() */
  118. /*---------------------------------------------------------------------------*/
  119. double
  120. _unur_dpdf_F(double x, const UNUR_DISTR *distr)
  121. {
  122. register const double *params = DISTR.params;
  123. if (x < 0.)
  124. /* out of support */
  125. return 0.;
  126. else if (_unur_iszero(x)) {
  127. if (nua < 2.)
  128. return -INFINITY;
  129. else if (_unur_isfsame(nub,2.))
  130. return -(2.+nub)/nub * exp(-LOGNORMCONSTANT);
  131. else
  132. return 0.;
  133. }
  134. else
  135. return _unur_pdf_F(x,distr) * _unur_dlogpdf_F(x,distr);
  136. } /* end of _unur_dpdf_F() */
  137. /*---------------------------------------------------------------------------*/
  138. double
  139. _unur_dlogpdf_F(double x, const UNUR_DISTR *distr)
  140. {
  141. register const double *params = DISTR.params;
  142. if (x < 0.)
  143. /* out of support */
  144. return 0.;
  145. else if (_unur_iszero(x)) {
  146. if (nua < 2.)
  147. return -INFINITY;
  148. else if (_unur_isfsame(nub,2.))
  149. return -(2.+nub)/nub;
  150. else
  151. return INFINITY;
  152. }
  153. else
  154. return ((nua/2.-1.)/x - nua*(nua+nub)/(2.*nub)/(1.+x*nua/nub));
  155. } /* end of _unur_dlogpdf_F() */
  156. /*---------------------------------------------------------------------------*/
  157. double
  158. _unur_cdf_F(double x, const UNUR_DISTR *distr)
  159. {
  160. #ifdef _unur_SF_cdf_F
  161. return _unur_SF_cdf_F(x,DISTR.nua,DISTR.nub);
  162. #else
  163. const double *params = DISTR.params;
  164. if (x <= 0.)
  165. /* out of support of p.d.f. */
  166. return 0.;
  167. if (nua * x > nub)
  168. return 1. - _unur_SF_incomplete_beta(nub / (nub + nua * x), nub/2., nua/2.);
  169. else
  170. return _unur_SF_incomplete_beta(nua * x / (nub + nua * x), nua/2., nub/2.);
  171. #endif
  172. } /* end of _unur_cdf_chisquare() */
  173. /*---------------------------------------------------------------------------*/
  174. #ifdef _unur_SF_invcdf_F
  175. double
  176. _unur_invcdf_F(double x, const UNUR_DISTR *distr)
  177. {
  178. return _unur_SF_invcdf_F(x,DISTR.nua,DISTR.nub);
  179. } /* end of _unur_invcdf_F() */
  180. #endif
  181. /*---------------------------------------------------------------------------*/
  182. int
  183. _unur_upd_mode_F( UNUR_DISTR *distr )
  184. {
  185. if (DISTR.nua >= 2.)
  186. DISTR.mode = ((DISTR.nua - 2.) * DISTR.nub) / (DISTR.nua * (DISTR.nub + 2.));
  187. else
  188. DISTR.mode = 0.;
  189. /* mode must be in domain */
  190. if (DISTR.mode < DISTR.domain[0])
  191. DISTR.mode = DISTR.domain[0];
  192. else if (DISTR.mode > DISTR.domain[1])
  193. DISTR.mode = DISTR.domain[1];
  194. return UNUR_SUCCESS;
  195. } /* end of _unur_upd_mode_F() */
  196. /*---------------------------------------------------------------------------*/
  197. int
  198. _unur_upd_area_F( UNUR_DISTR *distr )
  199. {
  200. /* log of normalization constant */
  201. LOGNORMCONSTANT = _unur_lognormconstant_F(DISTR.params,DISTR.n_params);
  202. if (distr->set & UNUR_DISTR_SET_STDDOMAIN) {
  203. DISTR.area = 1.;
  204. return UNUR_SUCCESS;
  205. }
  206. /* else */
  207. DISTR.area = ( _unur_cdf_F( DISTR.domain[1],distr)
  208. - _unur_cdf_F( DISTR.domain[0],distr) );
  209. return UNUR_SUCCESS;
  210. } /* end of _unur_upd_area_F() */
  211. /*---------------------------------------------------------------------------*/
  212. double
  213. _unur_lognormconstant_F(const double *params, int n_params ATTRIBUTE__UNUSED)
  214. {
  215. /* log( Beta(nu1/2, nu2/2) ) - (nu1/2) * log(nu1 / nu2) */
  216. return ((_unur_SF_ln_gamma(nua/2.) + _unur_SF_ln_gamma(nub/2.) - _unur_SF_ln_gamma((nua+nub)/2.))
  217. - 0.5 * nua * log(nua/nub));
  218. } /* end of _unur_lognormconstant_F() */
  219. /*---------------------------------------------------------------------------*/
  220. int
  221. _unur_set_params_F( UNUR_DISTR *distr, const double *params, int n_params )
  222. {
  223. /* check number of parameters for distribution */
  224. if (n_params < 2) {
  225. _unur_error(distr_name,UNUR_ERR_DISTR_NPARAMS,"too few"); return UNUR_ERR_DISTR_NPARAMS; }
  226. if (n_params > 2) {
  227. _unur_warning(distr_name,UNUR_ERR_DISTR_NPARAMS,"too many");
  228. n_params = 2; }
  229. CHECK_NULL(params,UNUR_ERR_NULL);
  230. /* check parameters nu1 and nu2 */
  231. if (nua <= 0. || nub <= 0.) {
  232. _unur_error(distr_name,UNUR_ERR_DISTR_DOMAIN,"nu <= 0");
  233. return UNUR_ERR_DISTR_DOMAIN;
  234. }
  235. /* copy parameters for standard form */
  236. DISTR.nua = nua;
  237. DISTR.nub = nub;
  238. /* copy optional parameters: none */
  239. /* store number of parameters */
  240. DISTR.n_params = n_params;
  241. /* set (standard) domain */
  242. if (distr->set & UNUR_DISTR_SET_STDDOMAIN) {
  243. DISTR.domain[0] = 0.; /* left boundary */
  244. DISTR.domain[1] = INFINITY; /* right boundary */
  245. }
  246. return UNUR_SUCCESS;
  247. } /* end of _unur_set_params_F() */
  248. /*---------------------------------------------------------------------------*/
  249. struct unur_distr *
  250. unur_distr_F( const double *params, int n_params )
  251. {
  252. register struct unur_distr *distr;
  253. /* get new (empty) distribution object */
  254. distr = unur_distr_cont_new();
  255. /* set distribution id */
  256. distr->id = UNUR_DISTR_F;
  257. /* name of distribution */
  258. distr->name = distr_name;
  259. /* how to get special generators */
  260. DISTR.init = NULL;
  261. /* functions */
  262. DISTR.pdf = _unur_pdf_F; /* pointer to PDF */
  263. DISTR.logpdf = _unur_logpdf_F; /* pointer to logPDF */
  264. DISTR.dpdf = _unur_dpdf_F; /* pointer to derivative of PDF */
  265. DISTR.dlogpdf = _unur_dlogpdf_F; /* pointer to derivative of logPDF */
  266. DISTR.cdf = _unur_cdf_F; /* pointer to CDF */
  267. #ifdef _unur_SF_invcdf_student
  268. DISTR.invcdf = _unur_invcdf_F; /* pointer to inverse CDF */
  269. #endif
  270. /* indicate which parameters are set */
  271. distr->set = ( UNUR_DISTR_SET_DOMAIN |
  272. UNUR_DISTR_SET_STDDOMAIN |
  273. UNUR_DISTR_SET_PDFAREA |
  274. UNUR_DISTR_SET_MODE );
  275. /* set parameters for distribution */
  276. if (_unur_set_params_F(distr,params,n_params)!=UNUR_SUCCESS) {
  277. free(distr);
  278. return NULL;
  279. }
  280. /* log of normalization constant */
  281. LOGNORMCONSTANT = _unur_lognormconstant_F(DISTR.params,DISTR.n_params);
  282. /* mode and area below p.d.f. */
  283. _unur_upd_mode_F( distr );
  284. DISTR.area = 1.;
  285. /* function for setting parameters and updating domain */
  286. DISTR.set_params = _unur_set_params_F;
  287. /* function for updating derived parameters */
  288. DISTR.upd_mode = _unur_upd_mode_F; /* funct for computing mode */
  289. DISTR.upd_area = _unur_upd_area_F; /* funct for computing area */
  290. /* return pointer to object */
  291. return distr;
  292. } /* end of unur_distr_F() */
  293. /*---------------------------------------------------------------------------*/
  294. #undef nu1
  295. #undef nu2
  296. #undef DISTR
  297. /*---------------------------------------------------------------------------*/