gravatar
alpar (Alpar Juttner)
alpar@cs.elte.hu
Pareto and Weibull random distributions
0 1 0
default
1 file changed with 28 insertions and 1 deletions:
↑ Collapse diff ↑
Show white space 16 line context
... ...
@@ -732,17 +732,17 @@
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) 
... ...
@@ -777,16 +777,43 @@
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
    
0 comments (0 inline)