Curve Fitting with Julia

Home

keywords: curve fitting, Weierstrass approximation, Hubbert peak theory, Seneca cliff, polynomial approximation

"How did you go bankrupt?” Bill asked.

"Two ways,” Mike said. “Gradually, then suddenly."

—Ernest Hemingway, The Sun Also Rises

The Weierstrass Approximation Theorem¶

In 1885, Karl Weierstrass

proved that any continuous function can be approximated as closely as you like with a polynomial. Written in mathematical terms, the theorem he proved is

Theorem: Suppose $f:[0,1]\Rightarrow \mathbb{R}$ is continuous. For any $ϵ>0$, there exists a polynomial $p$ such that

What does all this mean? Let's start with $f:[0,1]\Rightarrow \mathbb{R}$. This says that the function $f$ takes any value of $x$ between $0$ and $1$ and returns a real-valued number, $y=f(x)$ and $y$ can be anything. A function is continuous if it exists at every $x$ and you can trace the curve without lifting your pencil.

The part about "For any $ϵ>0$" is a trick mathematicians use a lot, and it means if you think of some very small positive number, we can always find another even smaller, but still positive number. The way we use this very, very tiny $ϵ$ is to take the absolute value of the difference between $f(x)$ and $p(x)$ at every $x$ and make sure it's less than $ϵ$.

The $\genfrac{}{}{0ex}{}{\text{sup}}{x\in [0,1]}$ is not "'Sup, dude?", but it's short for supremum and it means that for all possible polynomials $p$ approximating $f$, we want to choose the one where the maximum absolute value difference between $f(x)$ and $p(x)$ for any $x$ in the interval $[0,1]$ is the minimum. In other words, we want to choose the $p$ which minimizes the biggest error.

You might be wondering why the $x-$values are constrained to be between $0$ and $1$, and the answer is, it doesn't matter. You can shrink, stretch and shift $x$ to fit any interval. Suppose your function is defined over the interval $[a,b]$. Then, for every $x$, make

so when $x=a$, $x^{\prime }=0$ and when $x=b$, $x^{\prime }=1$. Then modify your function to

and you get the same $y$ values as $f$, but over the interval $[0,1]$.

A theorem isn't a theorem until it's been proved. Weierstrass proved his theorem, but in 1912, Sergei Bernstein

used Bernstein basis polynomials to construct approximation functions. In "Yacht Design with Mathematics" we used Bernstein polynomials as the basis for Bezier functions. For any $n\ge 0$, the $n+1$ Bernstein basis polynomials are defined as

with $x$ in the interval $[0,1]$. The Binomial Theorem shows how to expand the polynomial in $x$ and $y$ to the $n^{th}$ power,

Substitute $y=1-x$, add up all $k$ Bernstein polynomials and you get,

for any $x$. Bernstein then defined a new function which included a continuous function $f$ on the interval $[0,1]$ to be

Michael Bertrand gives the complete version of Bernstein's proof in his article, "Bernstein Proves the Weierstrass Approximation Theorem", and I'll provide another version at the end.

In 1938, Marshall Stone

extended Weierstrass' Theorem which is now known as the Stone-Weierstrass Theorem.

Some Applications¶

What can you do with a polynomial approximation to a function? With a good polynomial fit, you can use the polynomial to estimate the derivative or the integral of $f$, and this is what is done with software packages like Chebfun for Matlab, ApproxFun in Julia, or Python libraries such as pychebfun, ChebPy, or PaCAL. Polynomials are easy to differentiate and integrate, which is often not true for arbitrary functions, so having a good approximation is quite useful. These software tools can approximate continuous functions to machine precision, and computing with polynomials is fast.

With an approximate derivative or integral you can solve differential equations, and this has led to applications in neural networks as described in "Universal Invariant and Equivariant Graph Neural Networks". For insight into how Chebfun works, see "Chebfun: A New Kind of Numerical Computing" by Platte and Trefethen. The Chebfun website has many examples,

• Rational approximation of the Fermi-Dirac function describing the distribution of energy states for systems of particles.
• A drowning man and Snell's Law is derived from Fermat's principle giving the refraction of light through different media.
• Exponential, logistic, and Gompertz growth. "The Boonie Conundrum" used the Gompertz function to model the spread of a rumor.
• Procrustes shape analysis. We talked about the Greek myth of Procrustes in "Designing the Model".
• Modeling infectious disease outbreaks. In "How Deep Is That Pool?" we showed how to calculate the optimal pool size for COVID-19 testing.
• Smoothies: nowhere analytic functions. Weierstrass found functions that are continuous everywhere, but not differentiable anywhere.

You'll find plenty of nerdy-cool projects to try out. The only negative aspect of Chebfun is it requires the commercial software Matlab, but you may find an open-source version you prefer in Chebfun-related projects.

Polynomials aren't the only way to approximate functions. Fourier approximation and wavelets are other methods, and Nick Trefethen has written the book, "Approximation Theory and Approximation Practice" (partially available online). See the Publications page for other references.

Something you should never do is fit a function $f$ with a high degree polynomial $p$ and then extrapolate $p$ beyond the limits of $f$. You only get a good fit over the interval $[a,b]$ where you defined $f$.

Here, the function

may be approximated by the polynomial

Between $-\frac{\pi }{2}$ and $\frac{\pi }{2}$ the fit is pretty good, but outside the interval, the fit is much worse.

A Naive Approach to Approximations¶

Suppose you knew nothing about approximation theory and decided to try fitting one function (or data) with the sum of a simpler function?

Marion King Hubbert was a geologist, mathematician, and physicist who worked for Shell Oil, the USGS, and taught at Columbia, Stanford, and UC Berkeley. In 1956, he looked at oil production data and proposed a model for the amount of oil extracted from a field,

where $Q_{\text{max}}$ is the total quantity of oil to eventually be produced from the field, $Q(t)$ is the cumulative production up to time $t$, and constants $a$ and $b$ are parameters defined by local geology. This is the logistic function first described by Pierre Francois Verhulst to model population growth. The derivative of the logistic function produces a bell-shaped curve (or Gaussian function), which for Hubbert was the instantaneous rate of oil flow from the field.

Hubbert used the logistic function because initially with only a few wells production would be low, but as more came on line production should increase. Eventually, though, depletion would slow the rate of production. Hubbert thought his model could be extended to all U.S. oil production and predicted the peak of production would occur in 1970 for the lower 48 states.

In 1976 Hubbert predicted the world peak in production might happen as soon as 1995. A 2007 study by the UK Energy Resource Center found no geological reason to expect cumulative oil production should follow a logistic curve. They reported on field data suggesting peaks tended to happen at about 25% of the total production $Q_{\text{max}}$.

On the other hand, Ugo Bardi says decline is faster than growth and likes to quote Seneca,

"It would be some consolation for the feebleness of our selves and our works if all things should perish as slowly as they come into being; but as it is, increases are of sluggish growth, but the way to ruin is rapid." Lucius Anneaus Seneca, Letters to Lucilius, n. 91

We can model individual fields using Gaussian functions,

where $\alpha$ represents the peak amount produced, $\mu$ is the time of peak production and $\sigma$ determines the width of the bell curve. The Ghawar field in Saudi Arabia produces tens of thousands of barrels per day ($\alpha$), has been in operation since $\mu =1948$, and continues to this day so $\sigma$ is large. More recently, U.S. shale wells have increased production (the green curve departing from Hubbert's prediction), but their peaks are much lower and production drops 90% from the peak after only three years.

Is it possible to generate Bardi's Seneca Cliff with a series of Gaussians? Using Julia, the Seneca Cliff can be modeled as

f(x) = 0 < x <= 1 ? x : 0

which asks the question, Is $x$ between $0$ and $1$? If it is, $f(x)=x$, otherwise it's zero.

If oil production mimicked this $f(x)$, the rate would increase linearly until suddenly every oil well in the world stopped producing on the same day. This is about as likely as Zaphod Beeblebrox's Infinite Improbability Drive, but we're just trying to build a simple mind-sized problem. "See, all my procedures are mind-sized bites", explains seventh-grader Robert in Seymour Papert's Mindstorms.

First we need to create a vector of $x-$values,

x = Array(-0.5:0.01:1.5));

and then the corresponding $y-$values using a list comprehension

y = [f(x_i) for x_i in x];

Next, for each value of the parameters $\alpha ,\mu$ and $\sigma$ we can define a function $g$ and evaluate it over $x$,

g(x,α,μ,σ) = sum( α[i] * exp( - (x-μ[i])^2/(2*σ[i]^2) ) for i = 1:n)

or

which gives an estimate for $y$ using the sum of $n$ Gaussians

y_fit = [g(x_i,α,μ,σ) for x_i in x]

Using the LsqFit.jl package, we'll get the best parameters to fit $f(x)$. LsqFit requires an initial estimate, so I chose

and then created an initial vector $p_0=[\alpha ;\mu ;\sigma ];$ with $n=4$. To get the fit, LsqFit uses the function curve_fit

fit = curve_fit(gaussSum, x, y, p₀, lower = lb, upper = ub)

where gaussSum contains the definition for $g$, extracts the parameters from $p_0$, and generates $y\_\text{fit}$. The values of the parameters $\alpha ,\mu$ and $\sigma$ are required to be between the lower bound lb and the upper bound ub which were set to zero and one for this example.

The result is a pretty good approximation

You can see the peaks of each Gaussian with centers $\mu =[0.52,0.80,0.93,0.99]$ and widths $\sigma =[0.218,0.099,0.041,0.014]$. The optimal amplitudes are $\alpha =[0.51,0.57,0.63,0.69]$. The function $f$ has a discontinuity at $x=1$, so we wouldn't be able to find a polynomial approximation. Smoothing over the discontinuity would help with the sum of Gaussians technique. This doesn't mean this model necessarily represents what's actually happening, but it shows that a sum of Gaussian-like production curves can still lead to a Seneca Cliff. The Julia code is on GitHub.

This is a good experimental start to finding a method for fitting Gaussians to functions. Is there a pattern to the values for $\alpha ,\mu$ and $\sigma$? Can you find a functional representation for the parameters that works for different $n$? What happens if you change $f$? Can you find a general rule giving the parameters for any $f$ and any $n$? Test your hypotheses by changing $f$ (line 17) and plotting the fit.

The Proof¶

Here's the proof of Weierstrass' theorem. It's actually pretty easy to follow if you take it in mind-size bites. It comes from Chris Tisdell, a mathematician at UNSW Sydney. The proof below was taken from his YouTube video and is based on the proof by Dunham Jackson.

Introduction to the Theorem:

• Statement: Any continuous function defined on a closed interval $[a,b]$ can be uniformly approximated by a polynomial to any desired degree of accuracy.
• Importance: Establishes the foundation for polynomial approximation in analysis.

Key Concepts:

• Uniform Approximation: The concept that the maximum difference between the function and the approximating polynomial can be made arbitrarily small.
• Continuous Functions: Functions that are unbroken and have no gaps over the interval $[a,b].$

Constructing the Polynomial:

• Bernstein Polynomials: Introduction of Bernstein polynomials as a constructive method to approximate continuous functions.

• Definition: For a function $f(x)$ on $[0,1]$, the $n$-th Bernstein polynomial is given by:

  $$B_n(f,x)=\sum _{k=0}^{n}f\left(\frac{k}{n}\right)\left(\genfrac{}{}{0ex}{}{n}{k}\right)x^k(1-x)n-k.$$

Proof Outline:

• Step 1: Approximation on $[0,1]$:
  • Show that $B_n(f,x)$ converges uniformly to $f(x)$ on $[0,1]$ $n\to \infty$.
  • Use properties of binomial coefficients to demonstrate convergence for large $n$.
• Step 2: Extension to $[a,b]$:
  • Use a linear transformation to map the interval $[a,b]$ to $[0,1]$.
  • Apply the approximation on $[0,1]$ to the transformed function.
• Step 3: Uniform Convergence:
  • Prove that the maximum deviation $|f(x)-B_n(f,x)|$ can be made less than any $ϵ>0$ for sufficiently large $n$.
  • Use the properties of continuous functions and compactness of the interval to ensure uniform convergence.

Theorem: Suppose $f:[a,b]\Rightarrow \mathbb{R}$ is a continuous function. Then we can construct a sequence of polynomials $P_n$ such that for any $ϵ>0$ there exists $N=N(ϵ)$ such that for all $x\in [a,b]$ we have

Proof: Assume $0<a<b<1$. Extend $f$ to the entire real line by defining

This just adds pieces to $f(x)$ so that it covers the entire real line by connecting $(0,0)$ to $(a,f(a))$ and $(b,f(b))$ to $(1,0)$, and then defining $f(x)=0$ for any point outside the interval $[0,1]$, (and it's still continuous)

Now, create a sequence of constants (called a Landau integral),

Here are some examples of the Landau integral, and the numeric values over the interval $[-1,1]$. Integration gives the area below the curve and above the $x-$axis.

Next, define polynomials using these constants,

The first three polynomials of the sequence are

The idea is when $n$ gets to be big enough, $P_n(x)$ will very closely match $f(x)$. Notice that the degree of $P_n(x)$ is at most $2n$.

Since $f(x)=0$ for every $x$ outside of the interval $[0,1]$ then

Changing the limits of integration from 0 to 1 in the first line to $-1+x$ to $1+x$ in the second line doesn't really do anything for $x\in [0,1]$ except it now runs over $[-1,2]$ instead. In other words, let

so

where the first and third terms are zero since $f(x)$ is zero on those intervals, leaving the middle term unchanged.

If we change variables, letting $u=t-x$ so $t=u+x$, then $dt=du$. The limits of the integral change from

Now, look at the difference between this polynomial approximation $P_n(x)$ and the function we want to approximate $f(x)$,

where we use the Landau integral

Since $f$ is continuous, then for any $ϵ>0$ there is a $\delta$ such that

If $M$ is the maximum value of $|f|$ then

since $|f(x+u)|\le M$ and $|f(x)|\le M$, then

For any $|u|>\delta$, then $u^2>{\delta }^2$, so

and therefore

In either case, $|f(x+u)-f(x)|$ must be less than the sum of these two bounds, so

With these bounds, we can write

Let

and then using integration by parts ($\int \mu d\nu =\mu \nu -\int \nu d\mu$) and letting

So this means the integral $I_n$ can be bounded by $J_n/2n$,

In other words,

so long as $n$ is sufficiently large such that

Image credits¶

• Hero: Oil EIR Raises Sticky Issues For Hermosa Beach. Dana Murray, Heal the Bay, Apr. 14, 2014.
• Karl Weierstrass: Wikimedia Commons.
• Sergei Bernstein: Wikipedia.
• Marshall Stone: Marshall Harvey Stone seated in a classroom, ca. 1974. University of Massachusetts Amherst. Special Collections and University Archives, University of Massachusetts Amherst Libraries.
• Marion King Hubbert: Wikipedia.
• Hubbert Prediction Lower 48: Wikimedia Commons.
• The Seneca Cliff: The Seneca Effect. Ugo Bardi, The Seneca Effect: why decline is faster than growth. Aug 28, 2011.
• Dunham Jackson: Wikimedia Commons.

Code for this article¶

fitGaussParams.jl - Estimates parameters required to fit sum of Gaussians to a function: fitParams, gaussSum, plotFit

• Desmos - is an advanced graphing calculator implemented as a web application and a mobile application. In addition to graphing both equations and inequalities, it also features lists, plots, regressions, interactive variables, graph restriction, simultaneous graphing, piece wise function graphing, polar function graphing, two types of graphing grids – among other computational features commonly found in a programmable calculator.
• Julia - The Julia Project as a whole is about bringing usable, scalable technical computing to a greater audience: allowing scientists and researchers to use computation more rapidly and effectively; letting businesses do harder and more interesting analyses more easily and cheaply.
• VSCode and VSCodium - Visual Studio Code is a lightweight but powerful source code editor which runs on your desktop and is available for Windows, macOS and Linux. It comes with built-in support for JavaScript, TypeScript and Node.js and has a rich ecosystem of extensions for other languages and runtimes (such as C++, C#, Java, Python, PHP, Go, .NET). VSCodium is a community-driven, freely-licensed binary distribution of Microsoft’s editor VS Code.

References and further reading¶

• Bernstein Proves the Weierstrass Approximation Theorem, Michael Bertrand, Nonagon Math Library, Date Unknown.
• The Binomial Theorem, Wikipedia, Last accessed 04 Mar 2025.
• Chebfun: A New Kind of Numerical Computing, Platte, R. & Trefethen, L., SIAM Review, 2009.
• Universal Invariant and Equivariant Graph Neural Networks, Maron, H. et al., arXiv preprint, 2019.
• Approximation Theory and Approximation Practice, Trefethen, L.N., SIAM Books, 2013.
• Fourier Approximation, American Mathematical Society, 1997.
• The Logistic Function, Wikipedia, Last accessed 04 Mar 2025.
• The Seneca Effect: Why Decline is Faster Than Growth, Bardi, Ugo, Seneca Effect Blog, 28 Aug 2011.
• Shale Reality Check 2019, Hughes, David, Post Carbon Institute, 2019.
• Hubbert’s Peak: The Impending World Oil Shortage, Deffeyes, Kenneth, Princeton University Press, 2001.
• LsqFit.jl - Julia Nonlinear Least Squares Library, JuliaNLSolvers, GitHub Repository, Last accessed 04 Mar 2025.

📬 Subscribe and stay in the loop. · Email · GitHub · Forum · Facebook

© 2020 – 2025 Jan De Wilde & John Peach