gravatar
alpar (Alpar Juttner)
alpar@cs.elte.hu
Poisson distribution added
0 2 0
default
2 files changed with 24 insertions and 0 deletions:
↑ Collapse diff ↑
Ignore white space 6 line context
... ...
@@ -782,48 +782,71 @@
782 782
    /// Weibull distribution
783 783

	
784 784
    /// This function generates a Weibull distribution random number.
785 785
    /// 
786 786
    ///\param k shape parameter (<tt>k>0</tt>)
787 787
    ///\param lambda scale parameter (<tt>lambda>0</tt>)
788 788
    ///
789 789
    double weibull(double k,double lambda)
790 790
    {
791 791
      return lambda*pow(-std::log(1.0-real<double>()),1.0/k);
792 792
    }  
793 793
      
794 794
    /// Pareto distribution
795 795

	
796 796
    /// This function generates a Pareto distribution random number.
797 797
    /// 
798 798
    ///\param k shape parameter (<tt>k>0</tt>)
799 799
    ///\param x_min location parameter (<tt>x_min>0</tt>)
800 800
    ///
801 801
    double pareto(double k,double x_min)
802 802
    {
803 803
      return exponential(gamma(k,1.0/x_min));
804 804
    }  
805 805
      
806
    /// Poisson distribution
807

	
808
    /// This function generates a Poisson distribution random number with
809
    /// parameter \c lambda.
810
    /// 
811
    /// The probability mass function of this distribusion is
812
    /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f]
813
    /// \note The algorithm is taken from the book of Donald E. Knuth titled
814
    /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the
815
    /// return value.
816
    
817
    int poisson(double lambda)
818
    {
819
      const double l = std::exp(-lambda);
820
      int k=0;
821
      double p = 1.0;
822
      do {
823
	k++;
824
	p*=real<double>();
825
      } while (p>=l);
826
      return k-1;
827
    }  
828
      
806 829
    ///@}
807 830
    
808 831
    ///\name Two dimensional distributions
809 832
    ///
810 833

	
811 834
    ///@{
812 835
    
813 836
    /// Uniform distribution on the full unit circle
814 837

	
815 838
    /// Uniform distribution on the full unit circle.
816 839
    ///
817 840
    dim2::Point<double> disc() 
818 841
    {
819 842
      double V1,V2;
820 843
      do {
821 844
	V1=2*real<double>()-1;
822 845
	V2=2*real<double>()-1;
823 846
	
824 847
      } while(V1*V1+V2*V2>=1);
825 848
      return dim2::Point<double>(V1,V2);
826 849
    }
827 850
    /// A kind of two dimensional Gauss distribution
828 851

	
829 852
    /// This function provides a turning symmetric two-dimensional distribution.
Ignore white space 48 line context
... ...
@@ -12,25 +12,26 @@
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#include <lemon/random.h>
20 20
#include "test_tools.h"
21 21

	
22 22
///\file \brief Test cases for random.h
23 23
///
24 24
///\todo To be extended
25 25
///
26 26

	
27 27
int main()
28 28
{
29 29
  double a=lemon::rnd();
30 30
  check(a<1.0&&a>0.0,"This should be in [0,1)");
31 31
  a=lemon::rnd.gauss();
32 32
  a=lemon::rnd.gamma(3.45,0);
33 33
  a=lemon::rnd.gamma(4);
34 34
  //Does gamma work with integer k?
35 35
  a=lemon::rnd.gamma(4.0,0);
36
  a=lemon::rnd.poisson(.5);
36 37
}
0 comments (0 inline)