get the orthonormal Gegenbauer polynomials basis in R

I've been triying several packages to compute the orthonormal Gegenbauer polynomials basis which seems to be orthonormal at [-1,1]. Using the package midasmlI get however a different result:

library(midasml)

x <- seq(-1, 1, length.out = 100)
matplot(B <- gb(degree = 3, alpha = 1, a = -1, b = 1, jmax = NULL, X = x), type="l")

round(t(B) %*% B, 3) #not orthonormal

any reasonement why this is happening? any suggestions for an alternative approach? Thanks!

Solution

Orthonormality doesn't mean that this matrix is orthogonal. It means that the integral from -1 to 1 of P_i(x)P_j(x)w(x) is 0 if i != j and 1 if i = j, where w is a weight function. The weight function for the Gegenbauer polynomials with parameter alpha is (1 - x^2)^(alpha - 1/2).

You can get these polynomials with the orthopolynom package.

library(orthopolynom)
alpha <- 1
gp.list <- gegenbauer.polynomials(3, alpha, normalized = FALSE)
print(gp.list)
# [[1]]
# 1 
# 
# [[2]]
# 2*x 
# 
# [[3]]
# -1 + 4*x^2 
# 
# [[4]]
# -4*x + 8*x^3 

w <- function(x) { # the weight function
  (1 - x^2)^(alpha - 1/2)
}
GP2 <- function(x) {
  -1 + 4*x^2
}
GP3 <- function(x) {
  -4*x + 8*x^3
}

f <- function(x) {
  GP2(x) * GP3(x) * w(x)
}
integrate(f, lower = -1, upper = 1) # should be zero
# 0 with absolute error < 1.4e-14

f <- function(x) {
  GP3(x) * GP3(x) * w(x)
}
integrate(f, lower = -1, upper = 1) # should be pi/2
# 1.570796 with absolute error < 6e-07

Here the integral of P_i(x)²w(x) is not 1 because I took normalized = FALSE.