Best estimation of a fitting parameter to measured data

My goal is to estimate a parameter $$\alpha_1 = (\alpha_{11}, \alpha_{12})$$ which provides the best fit of certain measured data (a readout of some currents in a set of positions for a set of loads) to a hypothesized model.

After some hypothesis on the underlying data, the model to be fitted is represented by: $$I_{Ndc_{(1,k)}l} = \frac{I_{dc_{(k)}} - \alpha_{(1,k)} i_{ac_{(k)}}} {I_{dc_{(1)}} + I_{dc_{(2)}} - \alpha_{(1,1)} i_{ac_{(1)}} - \alpha_{(1,2)} i_{ac_{(2)}} }$$ where $$k=1,2$$. Therefore there are two functions, $$I_{Ndc_{(1,1)}l}$$ and $$I_{Ndc_{(1,2)}l}$$.

The measured data are the $$I_{dc_{(k)}}$$ and $$i_{ac_{(k)}}$$ and are represented by a $$n\times m$$ matrix, where each column is a different position and each row is a different load.

The idea is that there is a certain $$\alpha_1 = (\alpha_{11}, \alpha_{12})$$ which can make both $$I_{Ndc_{(1,k)}l}$$ independent of load variations (or independent to a variation of the 0 axis) for each position (or 1 axis).

Since I don't know how to formalize it and find the best estimator, I tried by intuition by finding (with python) the $$\alpha_1 = (\alpha_{11}, \alpha_{12})$$ which minimizes the unbiased variance (ddof=1) of $$I_{Ndc_{(1,k)}l}$$ in each position (means computing the variance on axis 0)

def res(alpha):
alpha_11, alpha_12 = alpha
INdc11l_meas = ((Idc1_meas - alpha_11 * iac1_meas)/
(Idc1_meas + Idc2_meas - alpha_11 * iac1_meas - alpha_12 * iac2_meas))

INdc12l_meas = ((Idc2_meas - alpha_12 * iac2_meas)/
(Idc1_meas + Idc2_meas - alpha_11 * iac1_meas - alpha_12 * iac2_meas))

var_INdc11l = sc.var(INdc11l_meas, axis=0, ddof=1) #Unbiased estimator ddof=1
var_INdc12l = sc.var(INdc12l_meas, axis=0, ddof=1)

return var_INdc11l + var_INdc12l

alpha_opt = optimize.minimize(res, [0,0])['x']


$$I_{Ndc_{(1,k)}l}$$ with alpha_opt results in superimposed lines (each line a different load):

What I am asking at this point is:

• Can this method be formalized in some way or have a particular meaning?
• Is there a better way to estimate $$\alpha_1$$?
• Is it possible to prove that it is the best estimator?
• What is the difference with minimizing the pairwise distance squared since it provides the same result? (changing the return statement with: return sc.sum(distance.pdist(INdc11l_meas, 'sqeuclidean')))

As a second step, with the $$\alpha_1$$ already known, it was needed to find a relation that can be inverted to find the position by measuring the currents $$I_{dc_{(k)}}$$ and $$i_{ac_{(k)}}$$. At this step they are not any more a matrix $$n \times m$$, but a scalar value, which provides a scalar $$I_{Ndc_{(1,k)}l}$$.

Again, by intuition, I supposed that by minimizing the variance, it is equivalent to minimize the distance of each point to the expectation and therefore the expectation to be a good candidate for that function

E = sc.mean(INdc11l_meas, axis=0)
f = sc.interpolate.interp1d(position_meas.T[:,0], E.T, 'cubic')
f_inv11 = lambda I: sc.optimize.root(lambda p,I: f(p) - I, [10], (I))

# From here Idc1_meas, Idc2_meas, iac1_meas, iac2_meas are scalars and instantaneous measurements
def find_pos11(Idc1_meas, Idc2_meas, iac1_meas, iac2_meas, alpha_1 = alpha_opt, f_inv11 = f_inv11):
return f_inv11(((Idc1_meas - alpha_11 * iac1_meas)/(Idc1_meas + Idc2_meas - alpha_1[0]* iac1_meas - alpha_1[1] * iac2_meas)))


The expectation is in yellow and the interpolated function in black:

What I am asking at this point is:

• Does it make sense?
• Is there a better way to estimate find_pos11?
• Is it possible to prove that it is the best available or to find which is the best?

Sorry in advance for my lack of rigour.

• I can provide a Python non-linear surface fitting example with both surface and contour plots if it might be of some use. Commented Oct 19, 2018 at 18:49
• Yes please. Which is the advantage over the method in the question? Commented Oct 19, 2018 at 18:59
• <<by intuition, I supposed that by minimizing the variance, it is equivalent to minimize the distance of each point to the expectation >> But how do you know where the expectation is? Commented Oct 19, 2018 at 19:02
• @user8948 The lines, for a certain $\alpha_1$, can be obtained using measured data and the expectation should be a line somewhere inside those lines, which ideally would be superimposed by hypotheses and therefore coincide with the expectation. Commented Oct 19, 2018 at 19:11
• Also: are your results good? That's important information. (I) "Will my results be good with this" and (II) "My results appear to be good, is this solid" are different questions. You can stress test your (appears-to-be-working) algorithm with cross-validation, the jackknife/bootstrap, etc. even if you can't math it out yet. Commented Oct 19, 2018 at 19:29

Per the comments, here is an example Python surface fitter using non-linear curve fitting with 3D scatter plot, 3D surface plot, and contour plot. This example uses manually supplied initial parameter estimates, but the scipy differential_evolution module has a genetic algorithm to help determine these initial estimates and I can give an example of that also.

import numpy, scipy, scipy.optimize
import matplotlib
from mpl_toolkits.mplot3d import  Axes3D
from matplotlib import cm # to colormap 3D surfaces from blue to red
import matplotlib.pyplot as plt

graphWidth = 800 # units are pixels
graphHeight = 600 # units are pixels

# 3D contour plot lines
numberOfContourLines = 16

def SurfacePlot(func, data, fittedParameters):
f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)

matplotlib.pyplot.grid(True)
axes = Axes3D(f)

x_data = data[0]
y_data = data[1]
z_data = data[2]

xModel = numpy.linspace(min(x_data), max(x_data), 20)
yModel = numpy.linspace(min(y_data), max(y_data), 20)
X, Y = numpy.meshgrid(xModel, yModel)

Z = func(numpy.array([X, Y]), *fittedParameters)

axes.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=1, antialiased=True)

axes.scatter(x_data, y_data, z_data) # show data along with plotted surface

axes.set_title('Surface Plot (click-drag with mouse)') # add a title for surface plot
axes.set_xlabel('X Data') # X axis data label
axes.set_ylabel('Y Data') # Y axis data label
axes.set_zlabel('Z Data') # Z axis data label

plt.show()
plt.close('all') # clean up after using pyplot or else thaere can be memory and process problems

def ContourPlot(func, data, fittedParameters):
f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)

x_data = data[0]
y_data = data[1]
z_data = data[2]

xModel = numpy.linspace(min(x_data), max(x_data), 20)
yModel = numpy.linspace(min(y_data), max(y_data), 20)
X, Y = numpy.meshgrid(xModel, yModel)

Z = func(numpy.array([X, Y]), *fittedParameters)

axes.plot(x_data, y_data, 'o')

axes.set_title('Contour Plot') # add a title for contour plot
axes.set_xlabel('X Data') # X axis data label
axes.set_ylabel('Y Data') # Y axis data label

CS = matplotlib.pyplot.contour(X, Y, Z, numberOfContourLines, colors='k')
matplotlib.pyplot.clabel(CS, inline=1, fontsize=10) # labels for contours

plt.show()
plt.close('all') # clean up after using pyplot or else thaere can be memory and process problems

def ScatterPlot(data):
f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)

matplotlib.pyplot.grid(True)
axes = Axes3D(f)
x_data = data[0]
y_data = data[1]
z_data = data[2]

axes.scatter(x_data, y_data, z_data)

axes.set_title('Scatter Plot (click-drag with mouse)')
axes.set_xlabel('X Data')
axes.set_ylabel('Y Data')
axes.set_zlabel('Z Data')

plt.show()
plt.close('all') # clean up after using pyplot or else thaere can be memory and process problems

def func(data, a, alpha, beta):
x = data[0]
y = data[1]
return a * (x**alpha) * (y**beta)

if __name__ == "__main__":
xData = numpy.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0])
yData = 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 = [xData, yData, zData]

# these are the same as scipy default values in this example
# if you do not have initial parameter estimates, they can
# be estimated using scipy's differential_evolution genetic
# algorithm module, I can supply an example if needed
initialParameters = [1.0, 1.0, 1.0]

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

ScatterPlot(data)
SurfacePlot(func, data, fittedParameters)
ContourPlot(func, data, fittedParameters)

print('fitted prameters', fittedParameters)

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)


It seems that minimizing the ensemble variance (as an objective function) is equivalent to minimizing the pair-wise distances between the variables

$$\text{FullSimplify}\left(\sum _{i=1}^N \left(y_i-\frac{\sum _{j=1}^N y_j}{N}\right)^2-\frac{\sum _{i=1}^N \sum _{j=i}^N \left(y_i-y_j\right){}^2}{N}\text{/.}\, \{N\to 20\}\right)=0$$

• Nice! I wonder how to extend to any N... Commented Oct 22, 2018 at 12:07