gravatar
alpar (Alpar Juttner)
alpar@cs.elte.hu
Serious buxfixes in Random::gamma() and Random::pareto()
0 1 0
default
1 file changed with 2 insertions and 2 deletions:
↑ Collapse diff ↑
Show white space 48 line context
... ...
@@ -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-gamma(int(std::floor(k))));
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)