blitz Version 0.9
|
00001 /*************************************************************************** 00002 * blitz/rand-normal.h Random Gaussian (Normal) generator 00003 * 00004 * $Id: rand-normal.h,v 1.4 2003/12/11 03:44:22 julianc Exp $ 00005 * 00006 * Copyright (C) 1997-2001 Todd Veldhuizen <tveldhui@oonumerics.org> 00007 * 00008 * This program is free software; you can redistribute it and/or 00009 * modify it under the terms of the GNU General Public License 00010 * as published by the Free Software Foundation; either version 2 00011 * of the License, or (at your option) any later version. 00012 * 00013 * This program is distributed in the hope that it will be useful, 00014 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00015 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00016 * GNU General Public License for more details. 00017 * 00018 * Suggestions: blitz-dev@oonumerics.org 00019 * Bugs: blitz-bugs@oonumerics.org 00020 * 00021 * For more information, please see the Blitz++ Home Page: 00022 * http://oonumerics.org/blitz/ 00023 * 00024 *************************************************************************** 00025 * 00026 * This generator transforms a (0,1] uniform distribution into 00027 * a Normal distribution. Let u,v be (0,1] random variables. Then 00028 * 00029 * x = sqrt(-2 ln v) cos(pi (2u-1)) 00030 * 00031 * is N(0,1) distributed. 00032 * 00033 * Reference: Athanasios Papoulis, "Probability, random variables, 00034 * and stochastic processes," McGraw-Hill : Toronto, 1991. 00035 * 00036 ***************************************************************************/ 00037 00038 #ifndef BZ_RAND_NORMAL_H 00039 #define BZ_RAND_NORMAL_H 00040 00041 #ifndef BZ_RANDOM_H 00042 #include <blitz/random.h> 00043 #endif 00044 00045 #ifndef BZ_RAND_UNIFORM_H 00046 #include <blitz/rand-uniform.h> 00047 #endif 00048 00049 #include <math.h> 00050 00051 BZ_NAMESPACE(blitz) 00052 00053 template<typename P_uniform BZ_TEMPLATE_DEFAULT(Uniform)> 00054 class Normal { 00055 00056 public: 00057 typedef double T_numtype; 00058 00059 Normal(double mean = 0.0, double variance = 1.0, double = 0.0) 00060 : mean_(mean), sigma_(::sqrt(variance)) 00061 { 00062 } 00063 00064 void randomize() 00065 { 00066 uniform_.randomize(); 00067 } 00068 00069 double random() 00070 { 00071 double u, v; 00072 00073 do { 00074 u = uniform_.random(); 00075 v = uniform_.random(); 00076 } while (v == 0); 00077 00078 return mean_ + sigma_ * ::sqrt(-2*::log(v)) * ::cos(M_PI * (2*u - 1)); 00079 } 00080 00081 private: 00082 double mean_, sigma_; 00083 P_uniform uniform_; 00084 }; 00085 00086 BZ_NAMESPACE_END 00087 00088 #endif // BZ_RAND_NORMAL_H 00089