How should I go about doing regression with a series of exponents involved? [duplicate]

Question

I have a table of data with input values and target values, and I was tasked to do something quite peculiar. I was tasked to run a sort of exponential regression and report back with the coefficients and the series of exponents. Let me explain:

Vanilla Linear Regression will minimize sum of squared error to create a linear combination in the format of $[a_1x_1 + a_2x_2...]$ where $a_i$ represents the series of coefficients and $x_i$ represents the different input features. This means that solving for the coefficients on linear regression is fairly easy given enough linearly independent data points.

I was given a dataset, and I was tasked to create an equation in the form of $[a_1x_1^{n_1} + a_2x_2^{n_2}...]$ where $a_i$ is the coefficients, $x_i$ is the input features, and $n_i$ signifies an exponential value.

I believe that since $a_i$ and $n_i$ need to be solved for, there are too many variables involved to go about this problem in a simple manner.

Is something like this even something like this even solvable? If so, how would you go about it statistically, and how would you go about it in python?

p.s. The task also says that all exponents, $n_i$, need to be between 0 and 1, so that might make it a bit easier.

p.s.s. I know that $[a_1x_1^{n_1} + a_2x_2^{n_2}...]$ is not linear unless all $n_i$ variables are either 0 or 1 but I wrote "linear regression" in the title because the two general formulas are similar.

Answer 1

Fitting to experimental data the function : $$y(x)=b\:x^{p}+c\:x^{q}$$ for the approximate values of the 4 parameters $p,q,b,c$ is easy thanks to the method described pp.71-72 (with numerical example p.73) in the paper https://fr.scribd.com/doc/14674814/Regressions-et-equations-integrales

The case of the 5 parameters function : $$y(x)=a+b\:x^{p}+c\:x^{q}$$ is shown in https://math.stackexchange.com/questions/3306953/sum-of-exponential-growth-and-decay/3317563#3317563

The extension to more parameters is theoretically possible on the same way. But too many parameters increase the deviations in the numerical calculus, making the result less and less accurate, up to be not acceptable in practice. This is discussed p.74 in the referenced paper.

Answer 2

Here is example code that I think should what you describe, using the equation "Z = (a1 * X1^n1) + (a2 * X2^n2)". This code uses the "brick wall" technique to keep the parameters n1 and n2 within bounds of 0.0 and 1 as noted in the comments in func().

fitted prameters:

a1 = 1.15465275488

n1 = 0.947251333112

a2 = 0.00698288550135

n2 = 0.999999491568

RMSE: 0.189049737197

R-squared: 0.99528439111

import numpy, scipy, scipy.optimize


def func(data, a1, n1, a2, n2):
    # if n1 or n2 is less than zero or greater than
    # one, return a large error - this will act as
    # a "brick wall" that the curve fit cannot overcome,
    # and this will act to keep the parameters within these
    # bounds. Note that the initial parameters estimates for
    # n1 and n2 (set below) must be within these bounds.
    if n1 < 0.0 or n1 > 1.0:
        return 1.0E300
    if n2 < 0.0 or n2 > 1.0:
        return 1.0E300

    x1 = data[0]
    x2 = data[1]
    return (a1 * x1**n1) + (a2 * x2**n2)


if __name__ == "__main__":
    x1Data = numpy.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0])
    x2Data = numpy.array([11.0, 12.1, 13.0, 14.1, 15.0, 16.1, 17.0, 18.1, 90.0])
    zData = numpy.array([1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.0, 9.9])

    data = [x1Data, x2Data, zData]

    # these are a1, n1, a2, n2 per the function definition
    # n1 and n2 are set within the bounds of 0.0 and 1.0
    initialParameters = [1.0, 0.5, 1.0, 0.5]

    # here a non-linear surface fit is made with scipy's curve_fit()
    fittedParameters, pcov = scipy.optimize.curve_fit(func, [x1Data, x2Data], zData, p0 = initialParameters)

    print('fitted prameters:')
    print('a1 =', fittedParameters[0])
    print('n1 =', fittedParameters[1])
    print('a2 =', fittedParameters[2])
    print('n2 =', fittedParameters[3])
    print()

    modelPredictions = func(data, *fittedParameters) 

    absError = modelPredictions - zData

    SE = numpy.square(absError) # squared errors
    MSE = numpy.mean(SE) # mean squared errors
    RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE
    Rsquared = 1.0 - (numpy.var(absError) / numpy.var(zData))
    print('RMSE:', RMSE)
    print('R-squared:', Rsquared)

Answer 3

If $x_i$ are the same

If you transform the $x$ values like:

$$t = \log(x)$$

Then your formula changes into

$$y = a_1 \exp(n_1t) + a_2 \exp(n_2t) + \dots$$

This you can solve with one of many methods for sums of exponential functions.

• An earlier example on Cross Validated is 'answer to: Fitting exponential decay with negative y values'.
• The article 'Exponential analysis in physical phenomena' by Andrei A. Istratov and Oleg F. Vyvenko published in Review of Scientific Instruments 70, 1233 (1999) provides a very thorough overview.

---

If $x_i$ are not the same

The above methods are for the case when the features $x_i$ are the same feature. E.g. a time variable. Your case might be a bit more difficult.

Solving $$y = a_1 \exp(n_1t) + a_2 \exp(n_2t) + \dots$$ is easier than $$y = a_1 \exp(n_1t_1) + a_2 \exp(n_2t_2) + \dots$$

For the latter, there is no easy way to find a good starting value. However, finding the solution can be made slightly easier by the method described below.

Segregated solving

You can solve it in a segregated fashion by considering that the function is a separable non-linear function and can be divided in a linear part and a non-linear part. When we try to minimize the loss function $\mathcal{L}(a_1,a_2,t_1,t_2)$ we can eliminate the coefficients $a_1$ and $a_2$ by considering them as the solution for an OLS fit with $t_1$ and $t_2$ fixed. Then these coefficients $a_1$ and $a_2$ are simple linear functions of $t_1$ and $t_2$ and we fill them into the equation for the loss function $\mathcal{L}(f_1(t_1,t_2),f_2(t_1,t_2),t_1,t_2)$ which is now effectively a function of only two parameters $\mathcal{L}(t_1,t_2)$. In this question an example is given for this method: How to compute the gradient for a seperable nonlinear least squares problem?

The segregated solving will optimize the algorithm to find the minimum, but you still need to find a starting point for $t_1$ and $t_2$. This is a difficult problem as the cost function is not convex. The example below (from the previously linked question) shows how the solution may diverge away from the global minimum.