Gamma function approximation

Question

I'm trying to write a three-parameter method that approximates the gamma function over a certain interval. The approximation should be a right-endpoint Riemann sum.

The gamma function is given by:

GAMMA(s) = 
inf
INT  x^(s-1) * exp(-x) dx
0

The right-endpoint Riemann sum approximation over the interval (0, m) should therefore be:

GAMMA(s) ~  
m
SUM  ((m/n)*i)^(s-1) * exp(-(m/n)*i) * delta_x        where delta_x = (m/n)
i=1

My code is as follows:

def gamma(x = 4.0, n = 100000, m = 2500)
  array = *(1..n)
  result = array.inject(0) {|sum, i| sum + ((((m/n)*i)**(x-1))*((2.7183)**(-(m/n)*i))*(m/n))}  
end

puts gamma

The code should return an approximation for 3! = 6, but instead it returns 0.0. Any ideas where I may be going wrong?

xlembouras · Accepted Answer · 2014-08-15T22:35:18.860

3

The problem is that when you do m/n

you are doing an integer division (eg 3/4 = 0) when you expect a float division (3/4 = 0.75)

you need to define your n and m as floats.

You can rewrite it as

def gamma(x = 4.0, n = 100000, m = 2500)
  n = n.to_f
  m = m.to_f

  (1..n).to_a.inject(0) do |sum, i|
    sum + ((((m/n)*i)**(x-1))*((Math::E)**(-(m/n)*i))*(m/n))
  end  
end

PS: also you do not need the array and result variables. PS2: consider using Math::E instead of 2.7183

edited Aug 15 '14 at 22:35

answered Aug 15 '14 at 22:29

xlembouras

8,215
4
33
42

I just implemented the process in excel and I got the correct answer, so I knew my approach was correct. I'm a coding n00b, so I figured it had to be that. Thanks again. This site is great. – Ice101781 Aug 15 '14 at 22:35

score 2 · Answer 2 · answered Aug 15 '14 at 23:24

2

Your problem was identified by @xlembouras. You might consider writing the method as follows.

Code

def gamma(x = 4.0, n = 100000, m = 2500)
  ratio = m.to_f/n
  xm1 = x-1.0
  ratio * (1..m).inject(0) do |sum,i|
    ixratio = i*ratio
    sum + ixratio**xm1 * Math.exp(-ixratio)
  end
end

Examples

gamma(x=4.0, n= 40, m =10).round(6) #=> 1.616233
gamma.round(6)                      #=> 6.0

Please confirm these calculations are correct.

answered Aug 15 '14 at 23:24

Cary Swoveland

106,649
6
63
100

gamma(4) should equal 3! = 6 – Ice101781 Aug 16 '14 at 17:22
I plan on writing a method that will extend the Gamma function to the entire complex plane. The ultimate goal is to be able to then write a method that will allow me to evaluate values of the functional equation of the Riemann Zeta function. – Ice101781 Aug 16 '14 at 17:25

Gamma function approximation

2 Answers2