fx-991 CW: Gamma Function

fx-991 CW: Gamma Function

All screenshots were made with Casio’s classpad.net website.

No Gamma Function? No Problem!

The gamma function is used a lot in advanced mathematics. The gamma functions with a lot of definition functions, but the most common one is for values t>0 (in particular Re(t)>0):

Γ(t) = ∫( x^(t-1) × e^(-x) dx, 0, ∞)

This integral is an improper integral and unless we have a calculator that handles infinite limit, we need to use the following:

Γ(t) =

lim ∫( x^(t-1) × e^(-x) dx, 0, w)

w → ∞

Calculators with the integral function can use the above for estimating the gamma function.

Use Basic Properties for Shortcuts

If t is a positive integer, we can use the factorial function:

Γ(t) = (t – 1)!

Example: Γ(16) = (16 – 1)! = 15! ≈ 1.31 × 10^12

If t is in the form n/2 where n is odd (i.e. 1/2 = 0.5, 3/2 = 1.5, 5/2 = 2.5, 7/2 = 3.5, etc.). we can use the product..

Γ(n / 2) = (n – 2) / 2 × (n – 4) / 2 × (n – 6) / 2 × … × 1 / 2 × √π

Example: Γ(3.5) = Γ(7 / 2) = 5 / 2 × 3 / 2 × 1 / 2 × √π = 15 / 8 × √π ≈ 3.32335097

The following pictures demonstrate the use of the above equivalency along with integral estimate:

I hope you find this helpful,

Eddie

 All original content copyright, © 2011-2025. Edward Shore. Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited. This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

The author does not use AI engines and never will.

HP 71B and HP Prime Python: Gaussian Quadrature

HP 71B and HP Prime Python: Gaussian Quadrature

A Very Brief Introduction

This is a top-level overview of the Gaussian Quadrature method. For more details, please check out the Sources section.

The design of GQ is to calculate:

∫( f(x) dx, -1, 1) ≈ Σ( wi * f(xi), i = 1, n)

Where do xi and wi come from?

* The xi values come from find the roots of the polynomials P’(x)= 0

* Pn’(x) are the derivatives of the Legendre polynomials Pn(x)

* After getting the roots, xi, the weights, w are calculated by the formula:

wi = 2 / ((1 -xi)^2 * Pn’(xi)^2)

Weights and Points of orders 3 and 5:

Among 3 points:

Roots (x) come from the derivative of the Legendre polynomial of the 3^rd order.

Pn’(x) = 3/2 * (5*x^2 – 1)

Point xi	Weight wi
-√(3/5) ≈ -0.7745 9669 2 = -√0.6	5/9
0	8/9
√(3/5) ≈ 0.7745 9669 2 = √0.6	5/9

Among 5 points:

Roots (x) come from the derivative of the Legendre polynomial of the 5^th order.

Pn’(x) = 15/8 * (21*x^4 – 14*x^2 + 1)

Point xi	Weight wi
-0.90617 98459	0.23692 68851
-0.53846 93101	0.47862 86705
0	128/225 ≈ 0.56888 88889
0.53846 93101	0.47862 86705
0.90617 98459	0.23692 68851

We can fit any interval [a, b] by this conversion:

∫( f(x) dx, a, b) = (b – a) / 2 * ∫( f( (b – a) / 2* x + (b + a) / 2 dx, -1, 1)

= (b – a) / 2 * Σ( wi * f( (b – a) / 2 * xi + (b + a) / 2 ), i = 1, n)

However, the programs on this blog entry will focus on the basic version of the Gaussian Quadrature.

The Code

The following code works with the following integrals:

∫( f(x) dx, -1, 1) ≈ Σ( wi * f(xi), i = 1, n)

HP 71B Basic

HP 71B: GQ3 (Gaussian Quadrature – 3 Point: Integrals from -1 to 1)

10 DEF FNF(X) = insert function of X here

15 S = 0

20 RADIANS

25 FOR I = 1 TO 3

30 READ X, W

31 DATA -SQR(0.6),5/9

32 DATA 0,8/9

33 DATA SQR(0.6),5/9

40 S = S + W * FNF(X)

55 NEXT I

50 DISP “ INTEGRAL= ”; S

HP 71B: GQ5 (Gaussian Quadrature – 5 Point: Integrals from -1 to 1)

10 DEF FNF(X) = insert function of X here

15 S = 0

20 RADIANS

25 FOR I = 1 TO 5

30 READ X, W

31 DATA -.9061798459, .2369268851

32 DATA -.5384693101, .4786286705

33 DATA 0, 128/255

34 DATA .5384693101, .4786286705

35 DATA .9061798459, .2369268851

40 S = S + W * FNF(X)

45 NEXT I

50 DISP “ INTEGRAL= ”; S

Integral Test: ∫( f(x) dx, -1, 1) ≈ Σ( wi * f(xi), i = 1, n)

Function f(x)	Actual (approx)	Gaussian 3 point (GQ3)	Gaussian 5 point (GQ5)
(x-3)*(x+4)*(x-5)	352/3 ≈ 117.333333333	117.333333333	117.333333334
-x^2 + 9*x + 10	58/3 ≈ 19.33333333	19.3333333334	19.3333333333
1.5 * cos(x)	2.524412954	2.52450532159	2.52441295561
exp(-x)	2.350402387	2.35033692868	2.35040238646
exp(sin(x))	2.283194521	2.28303911433	2.28319665353
sin(x)	0	0	0
1/(x - 2)	-1.09861228867	-1.09803921569	-1.09860924181

HP Prime Python App

HP Prime Python: gq3.py

# Gaussian Quadrature 3 point

from math import *

def f(x):

  # insert Function here

  return 1/(x-2)

s=0

x=[-sqrt(0.6),0,sqrt(0.6)]

w=[5/9,8/9,5/9]

for i in range(3):

  s=s+w[i]*f(x[i])

print("Integral: ",s)

HP Prime Python: gq5.py

# Gaussian Quadrature 5 point

from math import *

def f(x):

  # insert Function here

  return sin(x)

s=0

x=[0,-.5384693101056831,.5384693101056831,-.9061798459386640, .9061798459386640]

w=[.5688888888889,.4786286704993665,.4786286704993665,.2369268850561891,.2369268850561891]

for i in range(5):

  s=s+w[i]*f(x[i])

print("Integral: ",s)

Function f(x)	Actual (approx)	Gaussian 3 point (GQ3)	Gaussian 5 point (GQ5)
(x-3)*(x+4)*(x-5)	352/3 ≈ 117.333333333	117.333333333	117.333333333
-x^2 + 9*x + 10	58/3 ≈ 19.33333333	19.3333333333	19.3333333334
1.5 * cos(x)	2.524412954	2.52450532159038	2.52441295561081
exp(-x)	2.350402387	2.35033692868001	2.35040238646284
exp(sin(x))	2.283194521	2.2830391143268	2.28319665352905
sin(x)	0	0.0	0.0
1/(x - 2)	-1.09861228867	-1.09803921568627	-1.09860924181248

Sources

Kamermans, Mike “Pomax”. “Gaussian Quadrature Weights and Abscissae” https://pomax.github.io/bezierinfo/legendre-gauss.html 2011. Retrieved January 19, 2025.

Wikipedia. “Gauss-Legendre quadrature” https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrature Last edited January 19, 2025. Retrieved January 19, 2025.

Wikipedia. “Legendre Polynomials” https://en.wikipedia.org/wiki/Legendre_polynomials Last edited December 4, 2024. Retrieved January 19, 2025.

When it comes to evaluating integrals numerically, is Gaussian Quadrature better than Simpson’s Method?

Eddie

 All original content copyright, © 2011-2025. Edward Shore. Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited. This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

The content on this blog is 100% generated by humans. The author does not use AI engines and never will.

TI 30Xa Algorithm: The Sigmoid Function and its Integral

 TI 30Xa Algorithm: The Sigmoid Function and its Integral

The Sigmoid Function

The sigmoid function is defined as:

S = 1 / (1 + e^(-x))

If we multiply the last equation by e^(x) / e^(x), we get:

S = e^(x) / (e^(x) + 1) = e^(x) / (1 + e^(x))

The sigmoid function takes functions from the real numbers and maps them to the interval (0,1). If x → -∞, S → 0. If x → +∞, S → 1.

The inverse of the sigmoid function is the logit function, from which we can derive:

S = 1 / (1 + e^(-x))

1 / S = 1 + e^(-x)

1 / S – 1 = e^(-x)

ln( 1 / S – 1 ) = -x

ln( (1 – S) / S ) = -x (Note: 1 / S – 1 = 1 / S – S / S = (1 – S) / S )

-ln( (1 – S) / S ) = x

ln( S / (1 – S)) = x (Note: For any x, -ln(x) = ln(x^(-1)) = ln(1 / x))

The rest of the blog will focus on the sigmoid function.

The Integral of the Sigmoid Function

Funding the area of the curve under the sigmoid function is pretty straight forward.

∫ 1 / (1 + e^(-x)) dx

= ∫ e^(x) / (e^(x) + 1) dx

Let u = e^(x) + 1. The du = e^(x) dx and:

∫ du / u

= ln | u | + C

= ln (e^(x) +1) + C (since e^(x) + 1 > 0 for all real x)

The definite integral (area) can be calculated as:

x = b

∫ e^(x) / (e^(x) + 1) dx

x = a

= ln (1 + e^b) – ln (1 + e^a)

What if we are given the area, A (capital A), from negative infinity to a value x? Here we are finding the lower tail area.

t = x

∫ e^(t) / (e^(t) + 1) dx = A

t = -∞

ln(1 + e^(x)) – ln(1 + e^(-∞)) = A

Note that e^(-∞) → 0 as t → -∞

We can estimate that:

ln(1 + e^(x)) – ln( 1 ) = A

ln(1 + e^(x)) – 0 = A

1 + e^(x) = e^(A)

e^(x) = e^(A) – 1

x = ln(e^(A) – 1)

TI-30Xa Algorithms

Sigmoid function:

S = 1 / (1 + e^(-x)) = e^(x) / (e^(x) + 1)

Algorithm:

[ ( ] x [ +/- ] [ 2^nd ] [ LN ] {e^x} [ + ] 1 [ ) ] [ 1/x ] [ = ]

Example:

S(1.5):

[ ( ] 1.5 [ +/- ] [ 2^nd ] [ LN ] {e^x} [ + ] 1 [ ) ] [ 1/x ] [ = ]

Result: 0.817574476

Integral of the sigmoid function:

x = b

∫ e^(x) / (e^(x) + 1) dx

x = a

Algorithm:

[ ( ] b [ 2^nd ] [ LN ] {e^x} [ + ] 1 [ ) ] [ LN ]

[ - ] [ ( ] a [ 2^nd ] [ LN ] {e^x} [ + ] 1 [ ) ] [ LN ] [ = ]

Example:

x = 3

∫ e^(x) / (e^(x) + 1) dx

x = 0

(a = 0, b = 3)

[ ( ] 3 [ 2^nd ] [ LN ] {e^x} [ + ] 1 [ ) ] [ LN ]

[ - ] [ ( ] 0 [ 2^nd ] [ LN ] {e^x} [ + ] 1 [ ) ] [ LN ] [ = ]

Result: 2.355440171

Find x given lower tail area:

x = ln(e^(A) – 1)

Algorithm:

[ ( ] A [ 2^nd ] [ LN ] {e^x} [ - ] 1 [ ) ] [ LN ] [ = ]

Example:

Area: A = 0.5

[ ( ] 0.5 [ 2^nd ] [ LN ] {e^x} [ - ] 1 [ ) ] [ LN ] [ = ]

Result: -0.4327521296

Until next time,

Eddie

All original content copyright, © 2011-2024. Edward Shore. Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited. This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

Logit and Sigmoid Functions and its Calculus

Logit and Sigmoid Functions and its Calculus

Definitions

The sigmoid function is defined as:

sigmoid(x) = 1 ÷ (1 + e^(-x))

The logit function is defined as:  

logit(p) = ln (p ÷ (1 - p))

For logit(p) to have a real number answer, 0 ≤ p < 1

Transform from the Sigmoid Function to the Logit Function

We can easily transform from the sigmoid function to the logit function.  

Let s = sigmoid(x). Then:

s = 1 ÷ (1 + e^(-x))

s * (1 + e^(-x)) = 1

s + s * e^(-x) = 1

s * e^(-x) = 1 - s

e^(-x) = (1 - s) ÷ s

e^x = s ÷ (1 - s)

x = ln(s ÷ (1 - s)) = logit(s)

To transform from the logit function to the sigmoid function, just go backwards.  

Sigmoid Function:  Derivative and Integral

Derivative

s = sigmoid(x)

s = 1 ÷ (1 + e^(-x))

Using the quotient rule of derivatives:

ds/dx = [(1 + e^(-x)) * 0 - 1 * -e^(-x)] ÷ (1 + e^(-x))^2

= -(-e^(-x)) ÷ (1 + e^(-x))^2

= -e^(-x) ÷ (1 + e^(-x))^2

Integral

s = sigmoid(x)

s = 1 ÷ (1 + e^(-x))

Multiply both sides by e^x ÷ e^x:

s * (e^x ÷ e^x) = (e^x ÷ e^x) * (1 ÷ (1 + e^(-x)))

s = e^x ÷ (e^x + 1)

Integral:

∫ e^x ÷ (e^x + 1) dx

Let u = e^x + 1.  Then du = e^x dx 

= ∫  du ÷ (u + 1) 

= ln (u + 1) + C

= ln (e^x + 1) + C

Summary:

d/dx sigmoid(x) = -e^(-x) ÷ (1 + e^(-x))^2

∫ sigmoid(x) dx = ln (e^x + 1) + C

Logit Function:  Derivative and Integral

Derivative

logit(p) = ln (p ÷ (1 - p))

L = ln (p ÷ (1 - p))

Derivative:

dL/dp =  [(1 - p) ÷ p] * d/dp ln (p ÷ (1 - p))

=  [(1 - p) ÷ p] * [(1 - p) * 1 - p * (-1)] ÷ [(1 - p)^2] 

=  [(1 - p) ÷ p] * [1 - p + p] ÷ [(1 - p)^2]

=  [(1 - p) ÷ p] * 1 ÷ (1 - p)^2

= 1 ÷ [p * (1 - p)]

Integral:

∫ ln (p ÷ (1 - p)) dp

By integration by parts:

u = ln (p ÷ (1 - p)) 

du = 1 ÷ [p * (1 - p)] dp

v = dp

v = p

Then:

∫u dv

= p * ln ( p ÷ (1 - p)) - ∫ p ÷ (1 - p) dp

= p * ln ( p ÷ (1 - p)) + ∫ -p ÷ (1 - p) dp

= p * ln ( p ÷ (1 - p)) + ln(1 - p) + C

In Summary:

d/dp logit(p) = 1 ÷ [p * (1 - p)]

∫ logit(p) dp = p * ln ( p ÷ (1 - p)) + ln(1 - p) + C

Eddie

All original content copyright, © 2011-2022.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author. 

Python - Lambda Week: Integration by Simpson's Rule

Python - Lambda Week: Integration by Simpson's Rule

Welcome to Python Week!  This we we're going to cover calculus and the keyword lambda.

Note:  All Python scripts presented this week were created using a TI-NSpire CX II CAS.   As of June 2022, the lambda keyword is available on all calculators (in the United States) that have Python.   If you are not sure, please check your calculator manual. 

Simpson's Rule

The Simpson's Rule estimates numeric integrals by:

∫( f(x) dx, x = a to b) ≈

(b - a) /(3 * n) * (f(a) + 4 * f1 + 2 * f2 + 4 * f3 + .... + 2 * f_n-2 + 4 * f_n-1 + f(b))

n must be an even number of partitions.  The more partitions, the higher the accuracy and the higher computation time.

integrallam.py:  Numeric Integer

from math import *

print("The math module is imported.")

print("Integra of f(x), 6 places")

f=eval("lambda x:"+input("f(x)? "))

# input parameters

a=eval(input("lower = "))

b=eval(input("upper = "))

n=int(input("even parts: "))

# checksafe, add 1 if n is odd

if n/2-int(n/2)==0:

  n=n+1

# integral calculus

s=f(a)+f(b)

w=1

# 1 to n-1

for i in range(1,n):

  w=f(a+i*(b-a)/n)

  s+=(2*w) if (i/2-int(i/2)==0) else (4*w)

s*=(b-a)/(3*n)

print("Integral: "+str(round(s,6)))

All original content copyright, © 2011-2022.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.