r11:ea5945b2da9c
Fri, 21 Dec 2007 02:32:24
[Not Reviewed]
0 comments (0 inline)
0
1
0
28
1
| ... | ... |
@@ -716,93 +716,120 @@ |
| 716 | 716 |
double V1,V2,S; |
| 717 | 717 |
do {
|
| 718 | 718 |
V1=2*real<double>()-1; |
| 719 | 719 |
V2=2*real<double>()-1; |
| 720 | 720 |
S=V1*V1+V2*V2; |
| 721 | 721 |
} while(S>=1); |
| 722 | 722 |
return std::sqrt(-2*std::log(S)/S)*V1; |
| 723 | 723 |
} |
| 724 | 724 |
/// Gauss distribution with given mean and standard deviation |
| 725 | 725 |
|
| 726 | 726 |
/// \sa gauss() |
| 727 | 727 |
/// |
| 728 | 728 |
double gauss(double mean,double std_dev) |
| 729 | 729 |
{
|
| 730 | 730 |
return gauss()*std_dev+mean; |
| 731 | 731 |
} |
| 732 | 732 |
|
| 733 | 733 |
/// Exponential distribution with given mean |
| 734 | 734 |
|
| 735 | 735 |
/// This function generates an exponential distribution random number |
| 736 | 736 |
/// with mean <tt>1/lambda</tt>. |
| 737 | 737 |
/// |
| 738 | 738 |
double exponential(double lambda=1.0) |
| 739 | 739 |
{
|
| 740 |
return -std::log(real<double>())/lambda; |
|
| 740 |
return -std::log(1.0-real<double>())/lambda; |
|
| 741 | 741 |
} |
| 742 | 742 |
|
| 743 | 743 |
/// Gamma distribution with given integer shape |
| 744 | 744 |
|
| 745 | 745 |
/// This function generates a gamma distribution random number. |
| 746 | 746 |
/// |
| 747 | 747 |
///\param k shape parameter (<tt>k>0</tt> integer) |
| 748 | 748 |
double gamma(int k) |
| 749 | 749 |
{
|
| 750 | 750 |
double s = 0; |
| 751 | 751 |
for(int i=0;i<k;i++) s-=std::log(1.0-real<double>()); |
| 752 | 752 |
return s; |
| 753 | 753 |
} |
| 754 | 754 |
|
| 755 | 755 |
/// Gamma distribution with given shape and scale parameter |
| 756 | 756 |
|
| 757 | 757 |
/// This function generates a gamma distribution random number. |
| 758 | 758 |
/// |
| 759 | 759 |
///\param k shape parameter (<tt>k>0</tt>) |
| 760 | 760 |
///\param theta scale parameter |
| 761 | 761 |
/// |
| 762 | 762 |
double gamma(double k,double theta=1.0) |
| 763 | 763 |
{
|
| 764 | 764 |
double xi,nu; |
| 765 | 765 |
const double delta = k-std::floor(k); |
| 766 | 766 |
const double v0=M_E/(M_E-delta); |
| 767 | 767 |
do {
|
| 768 | 768 |
double V0=1.0-real<double>(); |
| 769 | 769 |
double V1=1.0-real<double>(); |
| 770 | 770 |
double V2=1.0-real<double>(); |
| 771 | 771 |
if(V2<=v0) |
| 772 | 772 |
{
|
| 773 | 773 |
xi=std::pow(V1,1.0/delta); |
| 774 | 774 |
nu=V0*std::pow(xi,delta-1.0); |
| 775 | 775 |
} |
| 776 | 776 |
else |
| 777 | 777 |
{
|
| 778 | 778 |
xi=1.0-std::log(V1); |
| 779 | 779 |
nu=V0*std::exp(-xi); |
| 780 | 780 |
} |
| 781 | 781 |
} while(nu>std::pow(xi,delta-1.0)*std::exp(-xi)); |
| 782 | 782 |
return theta*(xi-gamma(int(std::floor(k)))); |
| 783 | 783 |
} |
| 784 | 784 |
|
| 785 |
/// Weibull distribution |
|
| 786 |
|
|
| 787 |
/// This function generates a Weibull distribution random number. |
|
| 788 |
/// |
|
| 789 |
///\param k shape parameter (<tt>k>0</tt>) |
|
| 790 |
///\param lambda scale parameter (<tt>lambda>0</tt>) |
|
| 791 |
/// |
|
| 792 |
double weibull(double k,double lambda) |
|
| 793 |
{
|
|
| 794 |
return lambda*pow(-std::log(1.0-real<double>()),1.0/k); |
|
| 795 |
} |
|
| 796 |
|
|
| 797 |
/// Pareto distribution |
|
| 798 |
|
|
| 799 |
/// This function generates a Pareto distribution random number. |
|
| 800 |
/// |
|
| 801 |
///\param x_min location parameter (<tt>x_min>0</tt>) |
|
| 802 |
///\param k shape parameter (<tt>k>0</tt>) |
|
| 803 |
/// |
|
| 804 |
///\warning This function used inverse transform sampling, therefore may |
|
| 805 |
///suffer from numerical unstability. |
|
| 806 |
/// |
|
| 807 |
///\todo Implement a numerically stable method |
|
| 808 |
double pareto(double x_min,double k) |
|
| 809 |
{
|
|
| 810 |
return x_min*pow(1.0-real<double>(),1.0/k); |
|
| 811 |
} |
|
| 785 | 812 |
|
| 786 | 813 |
///@} |
| 787 | 814 |
|
| 788 | 815 |
///\name Two dimensional distributions |
| 789 | 816 |
/// |
| 790 | 817 |
|
| 791 | 818 |
///@{
|
| 792 | 819 |
|
| 793 | 820 |
/// Uniform distribution on the full unit circle. |
| 794 | 821 |
dim2::Point<double> disc() |
| 795 | 822 |
{
|
| 796 | 823 |
double V1,V2; |
| 797 | 824 |
do {
|
| 798 | 825 |
V1=2*real<double>()-1; |
| 799 | 826 |
V2=2*real<double>()-1; |
| 800 | 827 |
|
| 801 | 828 |
} while(V1*V1+V2*V2>=1); |
| 802 | 829 |
return dim2::Point<double>(V1,V2); |
| 803 | 830 |
} |
| 804 | 831 |
/// A kind of two dimensional Gauss distribution |
| 805 | 832 |
|
| 806 | 833 |
/// This function provides a turning symmetric two-dimensional distribution. |
| 807 | 834 |
/// Both coordinates are of standard normal distribution, but they are not |
| 808 | 835 |
/// independent. |
0 comments (0 inline)