Fitting power model to data with zeros

I want to fit a power model to some data. The model is of the form y = ax^b. I can take one of two approaches: 1) I can take logs to linearize the model, and use linear regression to find values for the parameters, or 2) I can use nonlinear regression.

My problem is that the first approach fails when some of the data are zero, and the second approach is sensitive to initial guesses of the parameters.

Which approach is better?

(Update: edited the sample data)

Here are some sample data:

x = [0, 20, 413, 454, 497, 598, 692]

y = [0, 0.4, 0.8, 1.3, 1.7, 2.1, 2.5]

Initially, I asked this question with these sample data (below), which should have an intercept, y = ax^b + c. This interests me, as well, but my question is really about the above set of sample data.

x = [0, 20, 413, 454, 497, 598, 692]

y = [135.0, 135.4, 135.8, 136.3, 136.7, 137.1, 137.5]

• As @James Phillips's reply indicates, it's hard to see these data as consistent with a power function of the kind you mention, as the data don't really suggest going through $0, 0$. Constant $y$ given $x$ is much more plausible than that. – Nick Cox Nov 2 '18 at 15:31

One solution is to replace the zero values with very tiny numbers, such as 1.0E-10. Here is a graphing Python example that does this with the scipy differential_evolution genetic algorithm module used to determine initial parameter estimates. I have added an offset to your equation, as this shifts the fitted model up and down the Y axis without changing the shape of the curve - the resulting fit seems much better in my opinion.

import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
import warnings

xData = numpy.array([0.0, 20.0, 413.0, 454.0, 497.0, 598.0, 692.0])
yData = numpy.array([135.0, 135.4, 135.8, 136.3, 136.7, 137.1, 137.5])

# replace zero values with extremely tiny values
numpy.place(xData, xData == 0.0, 1.0E-10)

def func(x, a, b, offset): # curve fitting function
return a*numpy.power(x, b) + offset

# function for genetic algorithm to minimize (sum of squared error)
def sumOfSquaredError(parameterTuple):
warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
val = func(xData, *parameterTuple)
return numpy.sum((yData - val) ** 2.0)

def generate_Initial_Parameters():
# min and max used for bounds
maxX = max(xData)
minX = min(xData)
maxY = max(yData)
minY = min(yData)

parameterBounds = []
parameterBounds.append([minY, maxY]) # search bounds for a
parameterBounds.append([minX, maxX]) # search bounds for b
parameterBounds.append([0.0, maxY]) # search bounds for Offset

# "seed" the numpy random number generator for repeatable results
result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
return result.x

# by default, differential_evolution completes by calling curve_fit() using parameter bounds
geneticParameters = generate_Initial_Parameters()

# now call curve_fit without passing bounds from the genetic algorithm,
# just in case the best fit parameters are aoutside those bounds
fittedParameters, pcov = curve_fit(func, xData, yData, geneticParameters)
print('Fitted parameters:', fittedParameters)
print()

modelPredictions = func(xData, *fittedParameters)

absError = modelPredictions - yData

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(yData))

print()
print('RMSE:', RMSE)
print('R-squared:', Rsquared)

print()

##########################################################
# graphics output section
def ModelAndScatterPlot(graphWidth, graphHeight):
f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)

# first the raw data as a scatter plot
axes.plot(xData, yData,  'D')

# create data for the fitted equation plot
xModel = numpy.linspace(min(xData), max(xData))
yModel = func(xModel, *fittedParameters)

# now the model as a line plot
axes.plot(xModel, yModel)

axes.set_xlabel('X Data') # X axis data label
axes.set_ylabel('Y Data') # Y axis data label

plt.show()
plt.close('all') # clean up after using pyplot

graphWidth = 800
graphHeight = 600
ModelAndScatterPlot(graphWidth, graphHeight)

• Very tiny replacements for 0 became large negative outliers on a logarithmic scale. I am not familiar with the methods you use, but can only say that that is a dubious suggestion for more mainstream methods, e.g. log transformation followed by least squares. – Nick Cox Nov 2 '18 at 1:15
• The graph proves that the method works in this case and obviously does not yield dubious results. As you see,there is no reason to perform log transform solely for the purpose of using linear algebra least squares regression, the non-linear regression without log transform works well - that is the main reason for including the plot and source code in the first place. – James Phillips Nov 2 '18 at 15:13
• Here "proves" is a strong word, but some sensitivity analysis is possible. In many cases where I've fitted a power function working on a logarithmic scale is really a good idea anyway. You've addressed the implausibility of a power function by adding an offset: could you address that in your answer? (I am in strong agreement that it's a much better model for these data.) – Nick Cox Nov 2 '18 at 15:29
• A scatterplot of the data indicated to me that an offset might be needed, so I tried it and got better results. The power function itself seems fine once the data offset is accounted for. – James Phillips Nov 2 '18 at 15:48