1

I've got some data points (X/Y coordinates) that were apparently created using a certain formula that I want to reconstruct now. I've only got those points, and I can plot them (e.g. like in this).
I know how to measure maxima, minima, slopes etc, and could solve an equation to compute the values for parameters in a formula so that it closely approximates/interpolates my data points, but which formula should I use?

I somehow need to find a parametric formula that matches the shape of the (parts of the) graph(s). What are the steps to identify such shapes, and to what formula (parts) do they correspond?

A step function is rather easy to identify, but I'm having problems with curves. Are they polynomials, of what degree? Are they trigonometric functions, like a sinus, squared, cubed? Do they contain a fraction, an exponentiation, something else? Or is it even something completely different, like a bezier curve?

I'm not looking for a solution to a specific problem, but a generic guide, hoping that this is not too broad.

MickG
  • 8,645
Bergi
  • 143

2 Answers2

3

The best way to handle this kind of problem of modelling data is to use information theoretic approaches.

Kolmogorov complexity

The simplest and completely generic way is to find a (program,input) pair with the shortest total length over a fixed general purpose programming language that when run outputs the data exactly as it is. This is called Kolmogorov complexity. Of course this depends on the chosen programming language, and clearly if you first look at the data you can construct a language purposely tailored to it (see Code Golf SE for some egregious examples). That's not allowed; you must first choose the language once and for all, and use it for every subsequent problem that you want to measure the complexity of the data. Also, the language must be Turing-complete. Common languages such as C++14 or Java 7 are suitable. The lengths of the program and input string is always measured in bits.

Let $K_L(x)$ be the Kolmogorov complexity of string $x$ over Turing-complete programming language $L$.

Why would this work? It can be easily proven that, given Turing-complete languages $L,M$, there is some $c$ such that, for any string $x$ we have $|K_L(x)-K_M(x)| \le c$. Thus for any class of strings the Kolmogorov complexity will be asymptotically correct as the string length goes to $\infty$. For example, the Kolmogorov complexity of the string $x_n$ of exactly $n$ zeros would differ from $\log_2(n)$ by at most $c$, for some constant $c$ that is independent of $n$. So as $n \to \infty$ we have $K_L(x_n) \in θ(\log(n))$.

Since our choice of language does not affect asymptotic Kolmogorov complexity, it is an objective measure of complexity for sufficiently large samples, and using a common language like C/C++ or Java makes it impossible for over-fitting.

Suitability of Kolmogorov complexity

For exact data, Kolmogorov complexity captures precisely the amount of information contained within it (up to a constant as mentioned above). However, it can be proven that there is no Turing-complete language $L$ and program $P$ such that $P(x) = K_L(x)$ for every string $x$. So you cannot hope to determine Kolmogorov complexity in general, not to say find the (program,input) pair that witnesses the complexity!

What you can still do is to find an upper bound, which is exactly what compression algorithms are essentially doing. Given any string $x$, the total length of the string that $x$ is compressed to and the decompression algorithm (written in $L$) is an upper bound on the true Kolmogorov complexity $K_L(x)$.

Nevertheless, most scientific experiments do not produce exact data, so it is pointless to ask for a program that reproduces the measurements exactly, since they are already imprecise and inaccurate. But the same ideas still apply. We just need to cater for noise.

Kolmogorov complexity modulo noise

If the data is a list of pairs $(x_k,y_k)_{k\in\{1..n\}}$, such as in many empirical experiments, and we know that the noise is uniform across the list, we can model the data as $N(P(x_k))=y_k$ where program $P$ is the actual function we seek and $N$ is a random program that models a noisy channel that adds noise. $N$ is written in a Turing-complete language that is additionally allowed to call a function that returns a uniformly random bit (which is independent of previous call results).

Then we want to a pair $(P,N)$ that minimizes what I call the complexity modulo noise defined as $|(P,N)| - \frac{1}{n} \sum_{k=1}^n \log_2(\mathbb{P}(N(P(x_k))=y_k))$. Intuitively, this is the length of the description of the trend and the noise source plus the average information in the noise itself. This is such that if there is no meaningful trend mapping $x$ to $y$, then the complexity modulo noise would be just the number of bits required to specify the noise. On the other hand if there is a trend then it would be captured by a short $P$, and there would be an optimal $N$ since it must be noisy enough to cater for the discrepancies but not more than necessary. Also, $N$ would capture the actual noise distribution if the sample is large enough.

Again, it is not possible to deterministically find $P,N$, but we can find upper bounds. Furthermore, this provides a way to compare approximations of an unknown function given a noisy sample of points. If we have 2 polynomials, one with low degree and simple coefficients that fits close enough, and the other with very high degree and complicated coefficients but passing through all the points at the sample resolution, you can be sure that the low degree polynomial will have lower complexity modulo noise. Likewise, you can try as many approximations as you like and simply choose the one that has the lowest complexity modulo noise.

This also allows you to use your intuition to guide your attempts. For real-world measurements, they tend to suffer from Gaussian noise, so you can try $N$ being a normal random variable with various variances (the mean is already captured by $P$).

It should be impossible to over-fit using this method, which is the crucial property that you are looking for.

Community
  • 1
user21820
  • 57,693
  • 9
  • 98
  • 256
  • Thanks, this is a very interesting approach. Unfortunately, it doesn't seem to come down to a simple "curve classification check list", but rather either brute-forcing or using intuition. – Bergi Apr 10 '16 at 18:10
  • @Bergi: That is one of my main points, which is that there cannot be a deterministic algorithm to classify a curve. However, it is not just about intuition, because this gives you a way to quantify how good your fit is. The complexity modulo noise tells you precisely how many bits your data intrinsically has, so if you find a very good fit, you can objectively compare its complexity against say the precision of the data itself. [continued] – user21820 Apr 11 '16 at 05:18
  • [continued] However, this method has a drawback due to the choice of programming language, which causes a high overhead that requires high precision in the data to offset. Another way I thought of is to define complexity modulo noise to be $|(P,N)| - \frac{1}{\sqrt{n}} \sum_{k=1}^n \log_2(\mathbb{P}(N(P(x_k))=y_k))$, so that it is as if we have a upper bound on the confidence interval for the actual amount of information, since the central limit theorem suggests that the standard deviation of what we measure would be proportional to $\sqrt{n}$. But I can't properly justify this method. – user21820 Apr 11 '16 at 06:27
  • @Bergi: Though I think the method in my above comment will work well since it heavily penalizes noise that is unnecessary, but if it is intrinsic then transferring noise into $(P,N)$ will incur greater cost. But there's a somewhat philosophical issue of how we should handle/ignore noise in experimental data when finding a model to fit it, which may be why I can't completely justify anything beyond the standard Kolmogorov complexity. – user21820 Apr 11 '16 at 06:36
  • You say "it can be easily proven that for any two Turing complete languages L, M, there is a constant c such that for any string x, |K_L(x) - K_M(x)| <= c. But this is false. Let L be C++ and let M be a language just like C++ except every character is doubled. i.e. in M you say iinntt mmaaiinn(()){{ instead of int main(){. The difference in lengths between programs in L and M is not bounded by a constant. – causative Feb 05 '22 at 16:52
  • @causative: You are completely wrong. Let C be an interpreter for L written in M. Let P be an M-program that on input (Q,x) runs Q via C on x. Then for every string x we have KL = |(Q,x)| for some L-program Q and also KM ≤ |(P,(Q,x))| ≤ |P|+O(log(|P|))+|(Q,x)| (using the obvious pair-coding). – user21820 Feb 05 '22 at 20:55
  • It's not right to say KL = |(Q, x)|. Isn't |(Q, x)| > |x|? That means you would have KL > |x| which is obviously not always the case (since it indicates no compression). Notation issues aside, I think you're trying to say that you can write an interpreter for L within M, called P, and then if Q is a shortest program in L that outputs x, then the program consisting of P and a quoted form of Q, call it (P,Q), would be a program in M that outputs x. That's true, but the problem is, the quoted form of Q within M is twice as long as Q, so (P,Q) within M is more than twice as long as Q. – causative Feb 05 '22 at 22:08
  • @causative: Sorry, I accidentally reused the variable "x", but my comment was in principle correct: Let C be an interpreter for L written in M. Let P be an M-program that on input (Q,x) runs Q via C on x. Then for every string x we have KL = |(Q,y)| for some L-program Q and string y and so KM ≤ |(P,(Q,y))| ≤ |P|+O(log(|P|))+|(Q,y)| (using the obvious pair-coding). Like I said, you need to actually read and think properly through my definition of K-complexity in the post before you claim it is wrong. – user21820 Feb 06 '22 at 15:02
  • There's still really no need for y. Anyway, like I said before, |(P, (Q, y))| within M, is about |P| + 2 |Q| + 2 |y|. Because you can't encode Q within M without doubling the length of Q. – causative Feb 06 '22 at 18:54
  • @causative: You're wrong; there is a need for y precisely to avoid having the stupid phenomenon you insist it has. If you want to not have y, you need to use prefix-free K-complexity, which some textbooks do, but that's besides the point here, which is that there is nothing wrong with my post. (Maybe you're still not reading it; I explicitly defined K-complexity in my post in terms of a (program,input) pair, and the input does not suffer from restrictions to the language!!) – user21820 Feb 06 '22 at 19:31
  • OK, I see what you're saying. – causative Feb 06 '22 at 19:48
  • @causative: Good, I'm glad that's cleared up. – user21820 Feb 06 '22 at 19:49
2

I do not know if I understood you correctly, but even if you want to find from a finite sequence of points (that have coordinates with a finite number of decimal places) in the plane whether the curve that goes through those points is polynomial or some non-polynomial continuous function it seems to me that it cannot be done because of the Stone–Weierstrass theorem which, roughly, says that you can uniformly approximate continuous function with polynomials to any degree of accuracy.

If it is the case that you have some finite set of points and you want to draw/find some, let us say, polynomial that goes through those points, then you have interpolation methods, see this one for example.

Farewell
  • 5,008
  • 1
    Thanks, yes, I know that I can approximate it using polynomials of high enough degree. But I'm looking for ways to establish that the curve was likely created from a polynomial at all, or maybe something else? If it was created with a polynomial, I want to restrict it to "simple" ones - that is, with few terms (many coefficients = 0) or of low degree. – Bergi Apr 10 '16 at 10:54