# Fitting a polynomial of degree 10 to data
**URL:** https://discourse.julialang.org/t/fitting-a-polynomial-of-degree-10-to-data/70924
**Category:** General Usage
**Tags:** data, regression, polynomials
**Created:** [November 4, 2021, 2:16am UTC](https://discourse.julialang.org/t/fitting-a-polynomial-of-degree-10-to-data/70924 "2021-11-04T02:16:27Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@valperez](https://discourse.julialang.org/u/valperez)
#### Post date: [November 4, 2021, 2:16am UTC](https://discourse.julialang.org/t/fitting-a-polynomial-of-degree-10-to-data/70924/1 "2021-11-04T02:16:27Z")
Hi there π
Iβm trying to test Juliaβs accuracy so Iβm trying to fit a polynomial of degree 10 to my data. The data can be found in the following link: [https://www.itl.nist.gov/div898/strd/lls/data/LINKS/DATA/Filip.dat](https://www.itl.nist.gov/div898/strd/lls/data/LINKS/DATA/Filip.dat)
I want to find the coefficients that best fit the problem
y = \sum\_{j = 0}^{j= 10 } beta\_j \* x^j + epsilon
I uploaded and named my data filip and I have a function called βpolyβ that creates a polynomial a certain degree.
So far I have tried:
1. fit(LinearModel, @formula(y ~ 1 + poly(x,10)), filip)
2. fm = @formula(y ~ 1 + poly(x, 10)) ; filip\_model = glm(fm, filip, Normal())
3. fm = @formula(y ~ 1 + poly(x, 10)); filip\_model = lm(fm, filip)
The problem is that I have the answer to the problem jaja is in this link: [StRD Certified Values for Dataset Filip](https://www.itl.nist.gov/div898/strd/lls/data/LINKS/v-Filip.shtml)
However, I want to know **how well** julia gets the decimals in the coeficients.
Thank you in advance π
---
### Author:  [@gustaphe](https://discourse.julialang.org/u/gustaphe)
#### Post date: [November 4, 2021, 4:23am UTC](https://discourse.julialang.org/t/fitting-a-polynomial-of-degree-10-to-data/70924/2 "2021-11-04T04:23:11Z")
Those are not base functions, which package are you using?
---
### Author:  [@nilshg](https://discourse.julialang.org/u/nilshg)
#### Post date: [November 4, 2021, 12:16pm UTC](https://discourse.julialang.org/t/fitting-a-polynomial-of-degree-10-to-data/70924/3 "2021-11-04T12:16:38Z")
Iβm assuming OP is using GLM with the additional `poly` function described in the docs [here](https://juliastats.org/StatsModels.jl/stable/internals/#extending).
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
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:
```julia
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.
---
### Author:  [@pdeffebach](https://discourse.julialang.org/u/pdeffebach)
#### Post date: [November 4, 2021, 12:21pm UTC](https://discourse.julialang.org/t/fitting-a-polynomial-of-degree-10-to-data/70924/4 "2021-11-04T12:21:07Z")
GLMβs `cholesky` call has a somewhat lower tolerance than Rβs. We need to add an option so users can change tolerance themselves in the `lm` call.
---
### Author:  [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)
#### Post date: [November 4, 2021, 1:25pm UTC](https://discourse.julialang.org/t/fitting-a-polynomial-of-degree-10-to-data/70924/5 "2021-11-04T13:25:03Z")
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.
---
### Author:  [@nilshg](https://discourse.julialang.org/u/nilshg)
#### Post date: [November 4, 2021, 2:11pm UTC](https://discourse.julialang.org/t/fitting-a-polynomial-of-degree-10-to-data/70924/7 "2021-11-04T14:11:22Z")
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.
---
### Author:  [@j\_verzani](https://discourse.julialang.org/u/j_verzani)
#### Post date: [November 4, 2021, 4:24pm UTC](https://discourse.julialang.org/t/fitting-a-polynomial-of-degree-10-to-data/70924/8 "2021-11-04T16:24:29Z")
Fitting using `Polynomials` is easy to do. I have `x`,`y` as variables and `as` as the given coefficients in the link:
```julia
using Polynomiales
p = fit(x,y, 10)
pa = fit(ArnoldiFit, x, y, 10)
q = Polynomial(as)
```
Then
```julia
julia> p - q
Polynomial(-4.5882096856075805e-5 - 8.77541174304497e-5*x - 7.449668601111625e-5*x^2 - 3.696358521665388e-5*x^3 - 1.187106710176522e-5*x^4 - 2.57854703988869e-6*x^5 - 3.836783157851187e-7*x^6 - 3.8622812281730035e-8*x^7 - 2.5178570339789985e-9*x^8 - 9.601549980434165e-11*x^9 - 1.6269660925026692e-12*x^10)
```
which are within the bounds given, on first glance. The ArnoldiFit is a better match to the given coefficients though:
```julia
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
```
---
### Author:  [@valperez](https://discourse.julialang.org/u/valperez)
#### Post date: [November 4, 2021, 5:27pm UTC](https://discourse.julialang.org/t/fitting-a-polynomial-of-degree-10-to-data/70924/9 "2021-11-04T17:27:03Z")
This worked!
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 π
---
### Author:  [@valperez](https://discourse.julialang.org/u/valperez)
#### Post date: [November 4, 2021, 6:17pm UTC](https://discourse.julialang.org/t/fitting-a-polynomial-of-degree-10-to-data/70924/11 "2021-11-04T18:17:44Z")
But they are very similar coefficients. I still have to figure out the standar deviation of the estimations.
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [November 4, 2021, 6:27pm UTC](https://discourse.julialang.org/t/fitting-a-polynomial-of-degree-10-to-data/70924/12 "2021-11-04T18:27:35Z")
Sorry my test seems to be bad one, will delete post.
_ **NB:**
Leave a picture just to illustrate the problem solved by_ @j_verzani:
