How to calculate confidence bands for models with 2 or more independent variables with kapteyn.kmpfit

Question

I try to compute confidence bands for a model with two independent variables.

I use kmpfit from kapteyn package to optimize parameters and compute confidence bands. The confidence_band function does not seem to work for models with more than one independent parameters. Does anyone have an idea how to solve this problem?

For models with one indpendent parameter, everything works as expected, without any problems:

from kapteyn import kmpfit
import numpy as np
import matplotlib.pyplot as plt 


#linear 1d model
def model_1d(p, x):
    a = p[0]
    b = p[1]
    return a*x + b


#data
x = np.array([0, 10, 20, 30])
y  = np.array([1.1, 1.8, 3.2, 3.0])

#optimizing parameters a and b with kmpfit
p0 = [0, 0] #start values for parameters a and b
f = kmpfit.simplefit(model_1d, p0, x, y)    


#calculation of confidence band

derivatives_a = x #1st derivative of a is x
derivatives_b = np.array([1, 1, 1, 1]) #1st derivative of b is always 1
dfdp = [derivatives_a, derivatives_b]

yhat, upper, lower = f.confidence_band(x, dfdp, 0.95, model_1d)


#plots
plt.plot(x,yhat) #plot fitted model
plt.plot(x,upper,'--') #plot upper confidence band
plt.plot(x,lower,'--') #plot lower confidence band
plt.plot(x,y,'.r') #plot data points

The script works nicely. As a result, you see the fitted model with upper and lower confidence bands: plot with upper and lower confidence bands

Below, the code is slightly changed to a 2D model, i.e. with two independent variables. This code does not work any more:

from kapteyn import kmpfit
import numpy as np

#linear with two independent variables x[0] amd x[1]
def model_2d(p, x):
    a = p[0]
    b = p[1]
    return a*x[0] + b*x[1]


#data
x = np.array([[0, 10, 20, 30], [0, 10, 20, 30]])
y  = [1.1, 1.8, 3.2, 3.0]

#optimizing parameters a and b with kmpfit
p0 = [0, 0] #start values for parameters a and b
f = kmpfit.simplefit(model_2d, p0, x, y)    


#calculation of confidence band

derivatives_a = x[0,:] #1st derivative of a is x[0]
derivatives_b = x[1,:] #1st derivative of b is x[1]
dfdp = [derivatives_a, derivatives_b]

yhat, upper, lower = f.confidence_band(x, dfdp, 0.95, model_2d)

For the 2D model, parameter optimization with simplefit still works. However calculation of confidence bands with confidence_band does not work any more. I get following error message, indicating that the shape of the numpy array x makes problems (however, the same x works for simplefit):

File "tmp.py", line 25, in <module>
    yhat, upper, lower = f.confidence_band(x, dfdp, 0.95, model_2d)
  File "src/kmpfit.pyx", line 845, in kmpfit.Fitter.confidence_band 
(src/kmpfit.c:8343)

If you have any idea how to solve this or if there are alternative python packages available, I would highly appreciate this. Thank you!

Answer 1

You might find lmfit (http://lmfit.github.io/lmfit-py/) useful for this. It has a robust but easy to use curve-fitting methods with its Model class, and also supports placing bounds and constraints on every Parameter.

In general, uncertainties for all variables and correlations between pairs of variables are calculated using the fast and generally reliable method of inverting the curvature matrix. For more detailed uncertainties (and including confidence of any specified 'sigma' level), it also allows for direct (if slightly slower) exploration of the confidence levels for any variable, and can help you make 2-d maps of chi-square for a pair of variables.

More specifically, I think what you are asking for is to generate the uncertainties in the model y from the uncertainties in the variables. Lmfit's Model class has an eval_uncertainty() method to do this. See http://lmfit.github.io/lmfit-py/model.html#calculating-uncertainties-in-the-model-function for an example.