1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
|
/*
$Id: scauchy.cc,v 1.9 2009/05/08 23:02:17 rhuey Exp $
AutoDock
Copyright (C) 2009 The Scripps Research Institute. All rights reserved.
AutoDock is a Trade Mark of The Scripps Research Institute.
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
/*
$-Id: scauchy.cc,v 3.0 1996/03/11 05:40:00 halliday Exp $
$-Source: /tmp_mnt/mgl/apps/src/autodock/3.0/autodock/RCS/scauchy.cc,v $
$-Log: scauchy.cc,v $
// Revision 3.0 1996/03/11 05:40:00 halliday
// The function definition for the GA/LS hybrid.
//
Revision 1.1 1993/02/05 15:42:11 whart
Initial revision
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
/* scauchy.cc
*
* Code for generating deviates from the Cauchy/Lorentzian distribution.
*
* Comments regarding this code are included at the end of the file.
*
* Code suggested by Lester Ingber <ingber@alumni.caltech.edu> and
* Numerical Recipies in C by Press et. al.
*/
#include <math.h>
#include "ranlib.h"
#define EPS 1.e-12
Real scauchy1()
{
Real x, y;
/* These four lines generate the tangent of a random
* angle; this is equivalent to
* y/x = tan(PI * ranf())
*/
do {
x = 2.0 * ranf() - 1.0;
y = 2.0 * ranf() - 1.0;
} while (x * x + y * y > 1.0);
if (fabs(x) < EPS)
x = (x < 0.0 ? x - EPS : x + EPS);
return (y / x);
}
Real scauchy2()
{
register Real x, y, r1, r2;
/* These four lines generate the tangent of a random
* angle; this is equivalent to
* y/x = tan(PI * ranf())
*/
do {
r1 = ranf();
// Using addition is faster than multiplication
x = r1 + r1 - 1.0;
r2 = ranf();
// Using addition is faster than multiplication
y = r2 + r2 - 1.0;
} while (x * x + y * y > 1.0);
if (fabs(x) < EPS)
x = (x < 0.0 ? x - EPS : x + EPS);
return (y / x);
}
/*******
Here's the summary of responses that I received concerning the
generation of Cauchy deviates:
Pretty much everyone who replied pointed out that the Cauchy distribution
can be generated by generating two standard normal deviates and taking
their quotient: N(0,1)/N(0,1). Barrett P. Eynon <barry@playfair.stanford.edu>
noted that the t distribution with 1 degree of freedom is equivalent to
the Cauchy distribution. Finally, B. Narasimhan <naras@cda.mrs.umn.edu>
noted that
the ratio X/Y is Cauchy any time the pair (X,Y) is radially
distributed. See Devroye, "Non Uniform Random Variate
Generation", for more. Thus you can even use ratio of two
coordinates of a point generated uniformly in the unit disk in
2D.
Code which uses this last observation was provided by Lester Ingber
<ingber@umiacs.UMD.EDU> and is included at the end of this article.
Mark Plutowski <pluto@cs.ucsd.edu> noted that the Cauchy distribution
is also called the Lorentzian distribution, and is discussed as such in
Numerical Recipies. Both NR and Ping Chou <choup@cgrb.orst.edu> noted
that you can generate a Cauchy deviate by taking a uniform deviate
on x \in [0,1] and computing tan(PI*x). However, in NR they later
compute the Cauchy deviates using the same method used in Ingber's code,
so presumably this is more efficient.
Robert E George <regeorge@magnus.acs.ohio-state.edu> also noted that
there is a routine in IMSL for generating Cauchy deviates.
Thanks for everyone's replies!
--Bill
*****/
|