I’m assuming OP is using GLM with the additional poly function described in the docs here.
Do I understand the question correctly as being about whether GLM.jl can replicate the coefficient values in given on the certification page you are linking? If so it appears it can’t:
julia> lm(@formula(y ~ x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + x^10), df)
StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}}}}, Matrix{Float64}}
y ~ 1 + x + :(x ^ 2) + :(x ^ 3) + :(x ^ 4) + :(x ^ 5) + :(x ^ 6) + :(x ^ 7) + :(x ^ 8) + :(x ^ 9) + :(x ^ 10)
Coefficients:
───────────────────────────────────────────────────────────────────────────────────────
Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95%
───────────────────────────────────────────────────────────────────────────────────────
(Intercept) 0.0 NaN NaN NaN NaN NaN
x 0.0 NaN NaN NaN NaN NaN
x ^ 2 0.0 NaN NaN NaN NaN NaN
x ^ 3 0.0 NaN NaN NaN NaN NaN
x ^ 4 0.0 NaN NaN NaN NaN NaN
x ^ 5 0.0 NaN NaN NaN NaN NaN
x ^ 6 0.00368644 0.000219028 16.83 <1e-26 0.0032503 0.00412258
x ^ 7 0.00191731 0.000127438 15.05 <1e-23 0.00166355 0.00217107
x ^ 8 0.000375885 2.75293e-5 13.65 <1e-21 0.000321067 0.000430703
x ^ 9 3.28191e-5 2.61667e-6 12.54 <1e-19 2.76087e-5 3.80296e-5
x ^ 10 1.07467e-6 9.23624e-8 11.64 <1e-17 8.90753e-7 1.25859e-6
───────────────────────────────────────────────────────────────────────────────────────
feeding the same data into R I get:
Call:
lm(formula = y ~ x + I(x^2) + I(x^3) + I(x^4) + I(x^5) + I(x^6) +
I(x^7) + I(x^8) + I(x^9) + I(x^10), data = df)
Residuals:
Min 1Q Median 3Q Max
-0.0099087 -0.0024610 0.0003385 0.0020743 0.0071654
Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.743e+02 8.756e+01 -1.990 0.050345 .
x -3.269e+02 1.480e+02 -2.208 0.030436 *
I(x^2) -2.661e+02 1.095e+02 -2.429 0.017617 *
I(x^3) -1.239e+02 4.652e+01 -2.664 0.009534 **
I(x^4) -3.638e+01 1.251e+01 -2.907 0.004845 **
I(x^5) -6.979e+00 2.211e+00 -3.156 0.002333 **
I(x^6) -8.747e-01 2.567e-01 -3.407 0.001079 **
I(x^7) -6.906e-02 1.890e-02 -3.654 0.000487 ***
I(x^8) -3.118e-03 8.009e-04 -3.894 0.000219 ***
I(x^9) -6.139e-05 1.489e-05 -4.123 9.91e-05 ***
I(x^10) NA NA NA NA
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.003768 on 72 degrees of freedom
Multiple R-squared: 0.9958, Adjusted R-squared: 0.9953
F-statistic: 1895 on 9 and 72 DF, p-value: < 2.2e-16
which seems somewhat closer but not quite there.
Maybe the problem is nearly rank-deficient, in which case GLM might behave differently from other software, but there might not necessarily be a “correct” approach to this issue always and everywhere. In this case it appears though that there is an “official” correct solution, maybe @nalimilan is interested in having a look at what’s happening here.
Fitting polynomials by solving for the coefficients of monomials is generally not a good idea. It generates poorly conditioned matrices and thus large errors in coefficients because the shapes of the curves x^n are very similar. Instead it’s better to express the polynomial as a linear combination of Chebyshev polynomials or Legendre polynomials, which are orthogonal to each other (in an inner product over the real interval and with proper weighting function, and at least strongly linearly independent when discretized).
This isn’t an esoteric problem. If you’re having getting good fits to a 10th order polynomial with order 10 data points, it’s likely your problem.
Not much point continuing the discussion without OP, but if I understand his post and the linked page correctly, this is a deliberately unstable problem which is used to validate statistical software against some known result. Now of course I don’t know how that known result has been derived and whether it’s in any way better than what R or GLM.jl give, but I think there’s little point arguing with the problem formulation.
which are within the bounds given, on first glance. The ArnoldiFit is a better match to the given coefficients though:
julia> maximum(abs(pa(x) - q(x)) for x ∈ range(0,1, length=1000))
7.275957614183426e-12
julia> maximum(abs(p(x) - q(x)) for x ∈ range(0,1, length=1000))
0.00025997101602115436
The only thing I changed was that I had to use p = Polynomials.fit(x,y, 10) and then q = Polynomial(as) and when I did the maximum(abs(p(x) - q(x)) for x ∈ range(0,1, length=1000)), the result gave me 0.0003723332774825394
I can’t tell you how grateful I am! thank you so so so much
