As regards extrapolation, I assume that this is an ill-posed inverse problem. However, it might well be that the maximum entropy approach can be applied to that problem. This is for instance used in recovering images that have been spoiled,e.g., by noise. [ http://www.et.byu.edu/~bjeffs/publications/Willis_ICIP_00.pdf ]

But this approach is not completely parameter free. Instead using a large parameter set, the bias is at least reduced by using that parameters that maximize the entropy. Of course, this does not help a lot if your model is too restricted even with the large parameter set.

I don't know too much about this. But the people in "nutonian" enterprise have a software for this http://www.nutonian.com/eureqa/research/

They describe their methods in this paper:

http://www.sciencemag.org/content/324/5923/81.short

Regards,

Nicolas

There is no "best" method. But if you are willing to assume a the function f(.) is reasonably smooth, then the cubic spline should give you a fairly good solution.

@Nicolas, I saw the relevant paper,and I think that it is a promising algorithmic process

@Lu, cubic splines interpolation can help only inside the given x-range of initial data, see attached Figure where outside [-1,1] fails to represent the true function.

I ask a method that can reproduce the true function AND outside the initial x-range.

Hi Demetris: If you want to estimate the function outside the range of observed xs, then it is essentially a extrapolation problem. Without additional knowledge or restriction on the function form, it is a non-identifiable problem from the statistical perspective. In any event, I think one may try cubic spline first to see the pattern of the function within the range of observed xs and then "guess" the reasonable extrapolation.

I suggest distance regression, compute the between units distance in X so to have all the N*(N-1)/2 distances in X you call Dx , then do the same thing for the same ordered couples in the Y space (by the way this procedure can be done for any number of X and Y variables). Now you can compute an ordinary linear regression between the Dx and Dy spaces, the obtaining of a strong correlation will tell you that a continuous transformation between the two spaces holds without imposing any special constraint different from continuity.....

It is not possible without having the 'domain' knowledge where the numbers are sampled from. In fact you can take it as a theorem and prove it easily.

Here I give a 'not-so-mathematical' proof - but you can formalize it easily from the basic argument.

Proof: Assume D={(x,y)} is the set of all x,y pair inputs you are supplying with-in the range. And C is the total corpus of pairs that is outside the range (D thus is essentially sampled from C).

Your goal is to predict C correctly from a method derived based on D.

Let us assume that you came-up with a mechanism or formalism F that observed D (and only D) and predicts C correctly.

The assumption is that F knows only D portion of C and does not know (or use) anything else about C.

Now, I take C and modify it slightly, but carefully, such that the D portion is not touched, but the post and pre parts are modified (perhaps by adding +1 to all y values for the same x values). Lets call this modified version as C1.

Both C and C1 are exactly same within the given input range of F - but out side that range they are completely different.

The question now is: Can your F predict C1 correctly if I just give the D portion of my C1??

Nope....unless it is a one-to-many function.

Thus, essentially you would need a non-deterministic automaton based functions to solve your problem - the non-determinism portion of which is determined by "extra domain knowledge" in the classical sense.

(Non-deterministic functions do not always return same output value for a given input. Their output vary each time you invoke the function, for the same input. But all our mathematics deal with deterministic functions with one-to-one repeat-ability)

Frankly if we can find one such F that can indeed predict out-of-range values from just the given range without extra knowledge, then all our problems will be solved in minutes. The question you are asking is much more bigger in scope than it appears and its solution can solve all problems, ranging from weather prediction, stock exchange, traffic behavior and so on.

Thank you,

Gopalakrishna Palem

http://blogs.msdn.com/gpalem

@Alessandro, after 'obtaining of a strong correlation will tell you that a continuous transformation between the two spaces holds' the problem still remains partially unsolved: We have not written out explicitly a functional form f(x) that can predict the data outside the xi \in [a,b]. I agree that we have established a qualitative kind of correlation and fruitful in information for the real yi=f(xi) function, but not the function itself.

@Gopalakrishna, your argument 'perhaps by adding +1 to all y values for the same x values' means a simple linear transformation and does not really change the functional form f(). I understand what you are trying to prove, but I have many counterexamples for your formal proof. Actually my work has shown me that 'as less is the inside information - such less is the accuracy of outside prediction' and not the non existence of such an option at all, as you claim.

Can you give us an example about ' a non-deterministic automaton based functions'?

The more you know, the better your chances are to find a reasonable model and/or a measure for your error. For instance, if you know that you have experimental data and you can derive a reasonable estimate of the error of each measurement, then your problem is better defined. Or if you know that all your data are positive and bounded everywhere while they vanish at infinity, then you might get even reasonable extrapolation curves if you are careful. If you don't know anything about the data or what "best" is in the context of the data, then you are in trouble. In my opinion, there is no unique a priori "best" method for the problem that you stated originally.

Demetris, in one sense what you are saying is true. The less you know about the data, the better the 'generalization'. I agree.

There is a catch, though: The less you know about the data, the more will be the 'cost' of the 'generalization'

There certainly are methods to work on 'generalized coordinates' - methods that consider the input-data with in the given range D, as 'generalized' coordinates (rather than as real data) and build the system for that 'generalized' data, that when transformed back to the real data solves real problems. But the cost of transforming to and from the 'generalized coordinates' is very costly.

In fact, some time back I was working on exactly such generalized methods for common problems (the attached paper gives some background) - but the main challenge is: the 'cost of transformation' to/from the generalized system out-weighs the regular solutions that use extra info about the data. (But still such methods will be useful for A.I. and Expert systems)

To come back to your argument: " 'perhaps by adding +1 to all y values for the same x values' means a simple linear transformation and does not really change the functional form f()"

That's correct: The new data C1 is just linear transformation of C (outside the range of D), with D being equal in both C and C1.

Now, I give you not one set C1, but many such data-sets C2, C3, C4 etc.... all share same portion of D, but different values out-side that range, by not just linear transformation, but some random polynomial transformation of some higher-degree.

The point is, if your fitting model F is correct, then it should just by observing D, correctly fit, not just one data set C, but *all* data sets C2, C3, C4 ......

The question is: how can F differentiate between C2 and C3 just by observing D (which is same in both). (We should not forget F, by definition is a static model and does not have any extra data information about any C_i).

In other words, you are trying to fit many functions into a single 'generic non-parametric function', the F.

In fact you can consider this F (the generic non-parametric method that can predict any given data-set) as the 'universal Turing Machine' mimicking the individual prediction functions f1, f2 .... fn, where each f_i is the correct prediction function for sequence C_i.

The 'parameter i' *is* the extra domain knowledge, that is required for such F to work.

Without the parameter that indicates which i should F mimic (being non-parametric), there is no way it can infer which state to go next from the given data D. That is the non-deterministic automaton.

That's where the contradiction is: "you are looking for a non-deterministic predictable function"

Non-deterministic methods, by definition, are un-predictable.

If you can predict something, then by definition, it becomes deterministic (because you have a mathematical procedure that tells exactly what the output is going to be for any given input).

As for working examples, there are no 'practical working methods' for non-deterministic automaton yet (again, because if we have a method, then it becomes deterministic !!!!)

But, you can understand it by assuming a function that has more than one value at any given point in its graph.

For example, a circle when you draw on a graph, it will have two y values for any single x value (except at, may be, extremes). So, given any x, the function will return + or - y. But which one to choose? That is the non-determinism. The method itself will not tell you which one to choose, its upto the 'viewer'.

If I ask someone to fit a model for such data (circular points) - which model fits best??

Sorry for the long reply - but the question you have raised is interesting one - in that it is on the converging border of many interesting theories from multiple fields - mathematics, computation theory, meta physics ...

Non-deterministic computation is very powerful - but unfortunately we are not, at least in our present day, harnessing any of its power in today's mathematics or computation - the prime reason being the underlying 'deterministic nature' of physics - the same reason why we cannot simulate 'true randomness' in computers.

Many open-challenges in the computational theory are just waiting for a good 'non-deterministic' method. I will not take more time - but just to give a brief of 'non-deterministic methods' - one of its applications is the so called 'quantum computing'

Understanding determinism/non-determinism is very crucial to many fields (as understanding the 'randomness' is crucial to probability) - even to statistics. After all, the Morkov Chains is indeed based on that (state-transformations)

You might like the 'determinism is ontic determinability is epistemic' paper by Harald Atmanspacher

Here is the my paper mentioned in the above answer.

Thank you,

Gopalakrishna Palem

http://gopalakrishna.palem.in

Article Data-dependencies and Learning in Artificial Systems

@Gopalakrishna, thanks for your extended reply. I think that I do not speak about the same 'task' (to speak in your familiar terminology). I do not want to find a function that should mimick all available instances of {(xi,yi),i=1,...,n} data set, but only to derive with the efficient time way and without adopting any 'model' the best underlying functional form that the specific data can give me. The quality of the extrapolation result is another story. Perhaps sometimes you cannot extrapolate due to singularities, for example for the function f(x)=1/(1+x^2) you cannot use its Taylor series expansion further of the interval (-1,1). So, my 'task' is to explore the available extrapolation techniques, those that they present good predictability.

The crucial point is the existence of extrapolation predictability and I have not solved this problem in general, because it is not possible to know every time something about the true nature of underlying data, as correctly @Herbert mentioned above.

As for the case of deterministic or not, we care about the situation where function f() EXISTS and we disregard other cases. Of course if some algorithm, like 'symbolic regression', mentioned from @Nicolas above, could judge if a 'law' simply does not exist, it could be a bonus, but we do not focus on it, we are searching for identifying the EXISTENT functional form.

As regards extrapolation, I assume that this is an ill-posed inverse problem. However, it might well be that the maximum entropy approach can be applied to that problem. This is for instance used in recovering images that have been spoiled,e.g., by noise. [ http://www.et.byu.edu/~bjeffs/publications/Willis_ICIP_00.pdf ]

But this approach is not completely parameter free. Instead using a large parameter set, the bias is at least reduced by using that parameters that maximize the entropy. Of course, this does not help a lot if your model is too restricted even with the large parameter set.

Dear Demetris, if the true function exists, I think the symbolic regression can be a good approach to discover it.

It is not easy to identify the correct function that originates the data (I agree with the comments form Gopalakrishna Palem) but by evolving several approximating functons you can capture the genotype of the true function. With genotype I mean the functional forms tha are persisting among all the evolved functions.

The following is my first attempt in this direction:

https://www.researchgate.net/publication/229084319_Understanding_Zooplankton_Long_Term_Variability_through_Genetic_Programming?ev=prf_pub

http://link.springer.com/chapter/10.1007%2F978-3-642-29066-4_5

I hope it helps

Simone

Conference Paper Understanding Zooplankton Long Term Variability through Gene...

@Simone, I think symbolic regression is the general way for solving this problem, although needs high computational resources. I have developed some other methods and I ' d like to compare my results with genetic programming results by using same data. (see my sent message for details).

As for the general problem, it is merely solved for the case of periodic data, f(x+T)=f(x), by using the proper Fourier analysis branch. But, the general, non periodic problem is still unsolved and probably it could not be well posed, as @Herbert noticed. I think that if that function really exists, then we have the ability to identify it in general, providing that we have reliable and sufficient data.

The problem is how can we define 'reliability' and 'sufficiency' for our data...

Gaussian processes are a classic such method, especially if you have some idea about the function's smoothness and noise.

More generally, a tree-based conditional density estimator ( such as one I developed a while back http://arxiv.org/abs/1005.2263 ) may be appropriate, if you have no idea what the noise may be.

I don't know there is a best way to do that. Any approache has its own limitation.

@Christos, I read your paper but I am not a 'Bayesian expert'. Is there any short worked example to present in brief the results of your method? Thank you.

Hi Demetris,

I suspect you're after a pragmatic answer to your question, so I'll suggest one. If your y is scalar (1-dimensional), or if it is n-dimensional but you feel itis safe to assume no interactions, then using generalised additive models (GAMs) is almost certainly close to ideal. If y multivariate (n-dimensional) and interactions among y's assumed, then random forests or similar (eg boosted regression trees; BRTs) or similar nearest neighbour based classifier /regression function ( e.g. Support vector machines; SVMs) or their variants are the current benchmark. Alternatively, if you feel it is safe to assume linear effects of y's on x, then some standard linear or generalised linear model (LM/GLM) are obviously sufficient.

Reputable and easy to use toolboxes are available in R for all of the above.

Of course, I'm speaking as an statistician. If you're interested in a answer, then I'll leave that to others to discuss.

Best wishes,

Matt

May be you can take a look to this article:

http://dx.doi.org/10.1111/j.1468-0394.2004.00278.x

"Evolving modular networks with genetic algorithms: application to nonlinear time series"

The idea implemented is the combination of a non-parametric regression model like the Artificial Neural Networks, with the knowledge of functional properties restrictions.

I hope it helps.

Regards

Antonio

You could take a look at Larry Wasserman's books All of Statistics: A Concise Course in Statistical Inference and All of Nonparametric Statistics (both Springer) to see if that helps (errata on LW's homepage at CMU Stat). Another extensive book is Sam Efromovich, Nonparametric Curve Estimation -- Methods, Theory, and Applications (Springer 1999). I agree that all this may be too much to go through, but it may give you further pointers. Best wishes -- Mihir

Mathew & Mihir, thank you for your suggestions. My main idea is to try solving such a problem not by the usual time series econometric apparatus but rather with other nonparametric methods, of course without excluding statistical concepts.

Antonio, the genetic algorithms are too promising, but computational expensive. Have you faced the problem of computational power in your research?(I am trying to find efficient methods, you can check my profile for details)

I think a light version of both methods is a suitable proposal.

Any other idea?

Demetris,

I have re-read you initial post. What you looking is non-parametric function regression. My answer .... a Multilayer Neural Network ..... but suffers for some overfitting if the parameters are no well established. Because training a neural networks is a non-linear optimization problem, you will face the same problems (non global minumun, .....)

This make the NN computationally expensive because you need to train many times, with different initial conditions and different subsets of original data.

Regards

Antonio

Hello,

Have you seen this article?

http://www.ssc.wisc.edu/~aaradill/zheng_1996.pdf

Hope this helps your research.

LPC

Unfortunately the too much promising genetic algorithms and symbolic regression techniques have been failed to recover a simple polynomial without errors imposed.

See for example:

>>>...

Total time elapsed: 26.552 seconds.

...

- G1DList

List size: 5

List: [20, -50, -50, -50, -2]

Source polynomial: [15, 14, 5, 12, -1]

Sample points: [(33, -748755), (7, 2073), (81, -36635475), (53, -6089155), (27, -291207), (21, -80835), (21, -80835), (-74, -34822905), (82, -38560977), (99, -84365607), (71, -21090535), (-93, -84415527), (1, 45), (19, -45927), (-97, -99435655), (14, -4297), (28, -346905), (59, -9634567), (62, -11896297), (39, -1593447)]

The initial polynomial was:

f(x) = 15+14*x+5*x^2+12*x^3-x^4

The recovered polynomial is:

g(x) = 20-50*x-50*x^2-50*x^3-2*x^4

They differ in an unacceptable way.

The test has been done with pyevolve package and python 2.6

http://pyevolve.sourceforge.net/ and http://www.python.org/

Any other and better solution in a reasonable time of execution?

@Demetris Christopoulos in your last answer you set the length of the chromosome to 5? so you already have an information about the structure of the function and you want to find the parameters ?

@Hossam, "I" am trying to recover a polynomial of 4th order, thus the reason why "I" have set the length equal to 5. Although this help, the algorithm over than 50% of the cases did not find even the correct sign for a[4]*x^4 term, leading to a polynomial towards the other side of the graph. Can you solve the problem with your setting and provide us the answer plus the time needed?

I am searching for the case of genetic algorithms, I am not an expert at that field.

Demetris, I have been using Random Forests for this type of problem with shockingly good results...see the book / papers by Breiman. Codes available in R and MATLAB.

This is fully non-parametric, but can also accommodate forests of "weak" learners that may have specified functional forms.

http://www.stat.berkeley.edu/~breiman/RandomForests/cc_home.htm

@Demetris I just used Cuckoo search which is a heuristic search algorithm to find the optimal coefficients of the 4th order polynomial depending on the 40 points which you used for GA.

The coefficients I got are [15.0956 13.9990 5.0001 12.0000 -1.0000 ] in exactly 30 seconds.

Check the the enclosed paper for details about Cuckoo Search algorithm

http://arxiv.org/pdf/1003.1594.pdf

@Hossam your result is excellent! Which implementation have you used for obtaining that result? (I found many based on Matlab, but not so user-friendly)

@John, can you suggest me an R package for Random Forests? I found many of those. Which is your preferred one?

@Demetris Please check the attached Matlab file which I customized to meet your polynomial and sample data. The file is ready to run in Matlab.

To change the Polynomial function => line 193

To change sample points => line 186-187

To change number of coefficients which you want to find => line 51

Demetris, I have been using MATLAB for this type of analysis, so I don't have a particular recommendation within R.

@Hossam, thank you very much!

@Hossam, I am trying to recover some data with cuckoo search, but it needs too much time for a problem that we could solve with a simple regression in seconds.

Try to find the 4th order polynomial with next data:

X= [-100, -95, -90, -85, -80, -75, -70, -65, -60, -55, -50, -45, -40, -35, -30, -25, -20, -15, -10, -5, 0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100];

Y=[-111951367.6, -91695295.62, -74318730.64, -59535163.58, -47073104.46, -36676047.90, -28102469.06, -21125909.64, -15534843.76, -11132737.43, -7738192.948, -5184623.179, -3320551.683, -2009476.317, -1129899.072, -575353.9909, -254251.3117, -90192.63870, -21610.83600, -2061.084832, 12.84106592, 1067.169579, 2642.084301, -8768.487678, -61711.76684, -199619.0606, -481080.2738, -979475.4633, -1783423.401, -2996346.723, -4736765.020, -7138223.486, -10349148.42, -14533076.41, -19868494.44, -26548922.27, -34782880.99, -44793807.88, -56820230.61, -71115672.73, -87948584.12];

I quit!

@Demetris, In your case we have more points so we need wider error margin. Try to change the Tolerance variable in line 51 to a larger value. Try to make it 400 or 360.

For me, I got [14.2704 14.1654 4.9984 12.0000 -1.0000] in less than 30 seconds for Tol=360.

Try it and update me please. Good luck !

P.S: Tolerance here is the summation of absolute errors.

@Hossam, thanks again:

>> tic;cuckoo_search_poly;toc

Total number of iterations=95400

fmin =

359.7510

bestnest =

15.5460 14.1558 4.9969 12.0000 -1.0000

Elapsed time is 6.196754 seconds.

(Tol=360)

This is great. Can you draw this polynomial compared to the optimal one please.

Yes, of course see Figure.

Some Remarks:

1)the method does not give all the times so close results

2)is there any formal way to choose the sum of absolute errors (tol) for every problem?

3)it is an alternative to regression, since you predefine the functional form and seek for the best parameters

Is there any generalization that can search alone for the best functional form, similar to symbolic regression?

*****************************

By the way: look a result of symbolic regression:

y = -x^4+12.0000436*x^3+5.00209460*x^2+108.0011772*x-2.93662250

It failed to recover the simple polynomial:

f(x) = -x^4+12*x^3+5*x^2+14*x+15

after the 3rd term.

@Demetris

regarding your 2nd point, usually there are two termination conditions of such evolutionary algorithm: 1) reaching minimum error value (tolerance) such as the one in our case 2) number of iterations is reached.

according to my knowledge, there is another cuckoo search algorithm which terminates after certain number of cycles. I will try to find it and configure it for you polynomial.

regarding your 3rd point, Genetic Programming (GP) and its variations are symbolic regression approaches which search for best functional form. You specify a set of simple mathematical operations like {*,-,+,/,sin,cos,log,power} and set of terminals (input variables) in our case is one variable X. Then GP tries to evolve a function using these sets that meet a set of given points while minimizing the error.

from my experience, its likely that GP will not retrieve your exact function but it can evolve one that gives good results according to the sample points (training points).

for example a GP run using same 41 sample points, GP retrieved the attached function with mean absolute error = 110.71.

Here is the retrieved function also in the image attached

@Hossam, I have used DTREG & pyevolve for GP recovery. The success is up to the user, if he/she can guess well, then many times the answer is good. But it is far away from concerning it as a 'non parametric method' (We' d say that it is a multi-parameter-searching method). Which is the software that you used for obtaining the above answer?

I use HeuristicLab framework. Download it, start the framework ( no installation needed). Then you can chose "Genetic Programming - Symbolic Regression" from the start page then load your sample points as a CSV file. Its easy to use.

Update me if you need further help.

http://dev.heuristiclab.com/trac/hl/core

OK, I used it and I cannot find a symbolic expression, just a solution tree, see Figure.

Is there any mandatory plugin which I have to install?

Thanks a lot.

OK, I found it!