... | ... |
@@ -692,141 +692,168 @@ |
692 | 692 |
bool boolean() { |
693 | 693 |
return bool_producer.convert(core); |
694 | 694 |
} |
695 | 695 |
|
696 | 696 |
///\name Nonuniform distributions |
697 | 697 |
/// |
698 | 698 |
|
699 | 699 |
///@{ |
700 | 700 |
|
701 | 701 |
/// \brief Returns a random bool |
702 | 702 |
/// |
703 | 703 |
/// It returns a random bool with given probability of true result |
704 | 704 |
bool boolean(double p) { |
705 | 705 |
return operator()() < p; |
706 | 706 |
} |
707 | 707 |
|
708 | 708 |
/// Standard Gauss distribution |
709 | 709 |
|
710 | 710 |
/// Standard Gauss distribution. |
711 | 711 |
/// \note The Cartesian form of the Box-Muller |
712 | 712 |
/// transformation is used to generate a random normal distribution. |
713 | 713 |
/// \todo Consider using the "ziggurat" method instead. |
714 | 714 |
double gauss() |
715 | 715 |
{ |
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. |
809 | 836 |
/// |
810 | 837 |
/// \note The coordinates are the two random variables provided by |
811 | 838 |
/// the Box-Muller method. |
812 | 839 |
dim2::Point<double> gauss2() |
813 | 840 |
{ |
814 | 841 |
double V1,V2,S; |
815 | 842 |
do { |
816 | 843 |
V1=2*real<double>()-1; |
817 | 844 |
V2=2*real<double>()-1; |
818 | 845 |
S=V1*V1+V2*V2; |
819 | 846 |
} while(S>=1); |
820 | 847 |
double W=std::sqrt(-2*std::log(S)/S); |
821 | 848 |
return dim2::Point<double>(W*V1,W*V2); |
822 | 849 |
} |
823 | 850 |
/// A kind of two dimensional exponential distribution |
824 | 851 |
|
825 | 852 |
/// This function provides a turning symmetric two-dimensional distribution. |
826 | 853 |
/// The x-coordinate is of conditionally exponential distribution |
827 | 854 |
/// with the condition that x is positive and y=0. If x is negative and |
828 | 855 |
/// y=0 then, -x is of exponential distribution. The same is true for the |
829 | 856 |
/// y-coordinate. |
830 | 857 |
dim2::Point<double> exponential2() |
831 | 858 |
{ |
832 | 859 |
double V1,V2,S; |
0 comments (0 inline)