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