Post Reply 
Handy Polynomial Fitting with Bernstein Polynomials
11-11-2018, 01:39 PM (This post was last modified: 11-12-2018 05:46 AM by Thomas Klemm.)
Post: #5
RE: Handy Polynomial Fitting with Bernstein Polynomials
I assumed you used the Bernstein Expansion to approximate a continuous function.

Here's a Python program that does this expansion:
Code:
from math import factorial

def binominal(n, k):
    return factorial(n) // factorial(k) // factorial(n-k)

def bernstein(n, k, x):
    return binominal(n, k) * x**k * (1-x)**(n-k)

def expand(f, x, n):
    return sum([f(j/n) * bernstein(n, j, x) for j in range(n+1)])

We can now compare the different results:
Code:
n = 10
for j in range(n+1):
    x = j / n
    f = lambda x: x**2
    approx = lambda x: x**2 + x*(1-x)/n
    print("%.1f: %.4f %.4f %.4f" % (x, f(x), expand(f, x, n), approx(x)))

x:   f(x)   expand approx
0.0: 0.0000 0.0000 0.0000
0.1: 0.0100 0.0190 0.0190
0.2: 0.0400 0.0560 0.0560
0.3: 0.0900 0.1110 0.1110
0.4: 0.1600 0.1840 0.1840
0.5: 0.2500 0.2750 0.2750
0.6: 0.3600 0.3840 0.3840
0.7: 0.4900 0.5110 0.5110
0.8: 0.6400 0.6560 0.6560
0.9: 0.8100 0.8190 0.8190
1.0: 1.0000 1.0000 1.0000


We can see that the expansion and the approximation agree.
That's astonishing since the expansion is based on Bernstein polynomials of degree 10 whereas the approximation is only of degree 2.

We can also see that approximation isn't really good except at \(x=0\) and \(x=1\).

However you appear to calculate a best fit of a function using Bernstein polynomials of a certain degree.

Here's how we can do this with Python:
Code:
from numpy import array, dot
from numpy.linalg import solve

b = array([x**2 for j in range(m+1) for x in [j/m]])
A = array([[bernstein(n, k, x) for k in range(n+1)] for j in range(m+1) for x in [j/m]])
These are the intermediate results:

>>> b
array([0. , 0.01, 0.04, 0.09, 0.16, 0.25, 0.36, 0.49, 0.64, 0.81, 1. ])

>>> A
array([[1.  , 0.  , 0.  ],
       [0.81, 0.18, 0.01],
       [0.64, 0.32, 0.04],
       [0.49, 0.42, 0.09],
       [0.36, 0.48, 0.16],
       [0.25, 0.5 , 0.25],
       [0.16, 0.48, 0.36],
       [0.09, 0.42, 0.49],
       [0.04, 0.32, 0.64],
       [0.01, 0.18, 0.81],
       [0.  , 0.  , 1.  ]])


For a best fit we solve the system that we get when both A and b are multiplied by A.T, the transpose of A:
Code:
solve(A.T @ A, A.T @ b)

The result is:

array([0., 0., 1.])

That's not a surprise as \(B_{2,2}(x) = x^2\).
However we would get the same result with fewer or more points as long as they are exactly on the curve as in Valentin's example.

What's the benefit of using Bernstein polynomials of order 2 compared to the ordinary basis \(1, x, x^2\)?

And why would you need an additional constant?
Since \(\sum_{i=0}^n B_{i,n}=1\) the constant 1 is already present as a linear combination.

Kind regards
Thomas
Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
RE: Handy Polynomial Fitting with Bernstein Polynomials - Thomas Klemm - 11-11-2018 01:39 PM



User(s) browsing this thread: 1 Guest(s)