... | ... |
@@ -774,73 +774,73 @@ |
774 | 774 |
/// |
775 | 775 |
///\param k shape parameter (<tt>k>0</tt>) |
776 | 776 |
///\param theta scale parameter |
777 | 777 |
/// |
778 | 778 |
double gamma(double k,double theta=1.0) |
779 | 779 |
{ |
780 | 780 |
double xi,nu; |
781 | 781 |
const double delta = k-std::floor(k); |
782 | 782 |
const double v0=E/(E-delta); |
783 | 783 |
do { |
784 | 784 |
double V0=1.0-real<double>(); |
785 | 785 |
double V1=1.0-real<double>(); |
786 | 786 |
double V2=1.0-real<double>(); |
787 | 787 |
if(V2<=v0) |
788 | 788 |
{ |
789 | 789 |
xi=std::pow(V1,1.0/delta); |
790 | 790 |
nu=V0*std::pow(xi,delta-1.0); |
791 | 791 |
} |
792 | 792 |
else |
793 | 793 |
{ |
794 | 794 |
xi=1.0-std::log(V1); |
795 | 795 |
nu=V0*std::exp(-xi); |
796 | 796 |
} |
797 | 797 |
} while(nu>std::pow(xi,delta-1.0)*std::exp(-xi)); |
798 |
return theta*(xi |
|
798 |
return theta*(xi+gamma(int(std::floor(k)))); |
|
799 | 799 |
} |
800 | 800 |
|
801 | 801 |
/// Weibull distribution |
802 | 802 |
|
803 | 803 |
/// This function generates a Weibull distribution random number. |
804 | 804 |
/// |
805 | 805 |
///\param k shape parameter (<tt>k>0</tt>) |
806 | 806 |
///\param lambda scale parameter (<tt>lambda>0</tt>) |
807 | 807 |
/// |
808 | 808 |
double weibull(double k,double lambda) |
809 | 809 |
{ |
810 | 810 |
return lambda*pow(-std::log(1.0-real<double>()),1.0/k); |
811 | 811 |
} |
812 | 812 |
|
813 | 813 |
/// Pareto distribution |
814 | 814 |
|
815 | 815 |
/// This function generates a Pareto distribution random number. |
816 | 816 |
/// |
817 | 817 |
///\param k shape parameter (<tt>k>0</tt>) |
818 | 818 |
///\param x_min location parameter (<tt>x_min>0</tt>) |
819 | 819 |
/// |
820 | 820 |
double pareto(double k,double x_min) |
821 | 821 |
{ |
822 |
return exponential(gamma(k,1.0/x_min)); |
|
822 |
return exponential(gamma(k,1.0/x_min))+x_min; |
|
823 | 823 |
} |
824 | 824 |
|
825 | 825 |
/// Poisson distribution |
826 | 826 |
|
827 | 827 |
/// This function generates a Poisson distribution random number with |
828 | 828 |
/// parameter \c lambda. |
829 | 829 |
/// |
830 | 830 |
/// The probability mass function of this distribusion is |
831 | 831 |
/// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f] |
832 | 832 |
/// \note The algorithm is taken from the book of Donald E. Knuth titled |
833 | 833 |
/// ''Seminumerical Algorithms'' (1969). Its running time is linear in the |
834 | 834 |
/// return value. |
835 | 835 |
|
836 | 836 |
int poisson(double lambda) |
837 | 837 |
{ |
838 | 838 |
const double l = std::exp(-lambda); |
839 | 839 |
int k=0; |
840 | 840 |
double p = 1.0; |
841 | 841 |
do { |
842 | 842 |
k++; |
843 | 843 |
p*=real<double>(); |
844 | 844 |
} while (p>=l); |
845 | 845 |
return k-1; |
846 | 846 |
} |
0 comments (0 inline)