Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Evaluating the Bell polynomials #139

Open
ignace-computing opened this issue Jun 29, 2023 · 1 comment
Open

Evaluating the Bell polynomials #139

ignace-computing opened this issue Jun 29, 2023 · 1 comment

Comments

@ignace-computing
Copy link

Dear contributors of the Combinatorics.jl package,

Would you please know where I could find an implementation of (the evaluation) of the incomplete Bell polynomials?
https://en.wikipedia.org/wiki/Bell_polynomials#Definitions#Exponential%20Bell%20polynomials

Thank you!

@ignace-computing
Copy link
Author

Hello,

I just wrote a function, based on a recurrence relation that can be found on Wikipedia.

Do you have any suggestion whether it is useful for this package, or else if this code could be a useful addition elsewhere?
Any (other) comment of course welcome too.

Cheers,

Evaluate the incomplete Bell polynomial B_{N,K}(x)
where x is an array containing the abscissa values [x_1, ... x{n-k+1}]

Definition and background information, see 
https://en.wikipedia.org/wiki/Bell_polynomials#Definitions
https://en.wikipedia.org/wiki/Bell_polynomials#Properties
https://en.wikipedia.org/wiki/Bell_polynomials#Table_of_values
"""
evaluate_incomplete_Bell_poly = function(N,K,all_x)
	# see formula C2
	# let us construct the triangular table
	if K > N
		return 0.
	end
	all_B_old = zeros(N+1,1)
	all_B_old[1] = 1.
	all_B_new = zeros(N+1,1)
	all_B_new[1] = 1.
	for k=1:K
		for n=1:N
			all_B_new[n+1] = sum([binomial(n-1,i-1)*all_x[i]*all_B_old[n-i+1] for i=1:n-k+1])
		end
		all_B_old .= all_B_new
	end
	return all_B_new[end]
end

## test the function 
x = randn(11,1)

@assert evaluate_incomplete_Bell_poly(0,0,x) ≈ 1.

@assert evaluate_incomplete_Bell_poly(1,0,x) ≈ 0.
@assert evaluate_incomplete_Bell_poly(3,0,x)≈ 0.

@assert evaluate_incomplete_Bell_poly(0,1,x) ≈ 0.
@assert evaluate_incomplete_Bell_poly(0,3,x) ≈ 0.

@assert evaluate_incomplete_Bell_poly(1,1,x) ≈ x[1]
@assert evaluate_incomplete_Bell_poly(2,2,x) ≈ x[1]^2
@assert evaluate_incomplete_Bell_poly(3,3,x) ≈ x[1]^3

@assert evaluate_incomplete_Bell_poly(6,4,x) ≈ 45*x[1]^2*x[2]^2 + 20*x[1]^3*x[3]
@assert evaluate_incomplete_Bell_poly(5,4,x) ≈ 10*x[1]^3*x[2]
@assert evaluate_incomplete_Bell_poly(6,3,x) ≈ 15*x[2]^3 + 60*x[1]*x[2]*x[3] + 15*x[1]^2*x[4]

# use some results in (Abbas, M.; Bouroubi, S. (2005). "On new identities for Bell's polynomial". Discrete Math. 293 (1–3)) as reference
# see https://vinar.vin.bg.ac.rs//bitstream/id/12791/4374.pdf
@assert evaluate_incomplete_Bell_poly(9,7,x) ≈ 378*x[1]^5*x[2]^2 + 84*x[1]^6*x[3]
@assert evaluate_incomplete_Bell_poly(10,7,x) ≈ 3150*x[1]^4*x[2]^3 + 2520*x[1]^5*x[2]*x[3] + 210*x[1]^6*x[4]
@assert evaluate_incomplete_Bell_poly(11,7,x) ≈ 17325*x[1]^3*x[2]^4 + 34650*x[1]^4*x[2]^2*x[3] + + 4620*x[1]^5*x[3]^2 + 6930*x[1]^5*x[2]*x[4] + 462*x[1]^6*x[5]

println("tests passed")

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant