user9259974
user9259974

Reputation:

How to fit data with two set of independent variables

i have this equation a*(t^alpha)*(p_p^beta) that i want to fit to get the alpha and beta values where t and p_p are the independent variables. My questions is how do i write the final fit model (result) expression.

result = model.fit(S_L1, params, t=t, p_p=p_p)

I tried something like the above expression but i am getting this error:

ValueError: The input contains nan values


# Calculating unburned mass temperature
T_u = T_i*(p_filter/p_i)**((k_u-1)/k_u)                     # Linear unburned temperature
t = T_u/T_i
p_p = p_filter/p_i


# Model function.
def mod_m(t, p_p, a=1, alpha=1,beta=1):                             # Define function with initial guesses
    return a*(t**alpha)*(p_p**beta)                                 # Function for fitting


# Fitting model.
model = Model(mod_m, independent_vars=['t','p_p'] )

# Making a set of parameters:
params = model.make_params(a=10)

# Setting  min/max bounds on parameters:
params['alpha'].min = 0.0
params['beta'].min = 0.0
params['a'].min = 0.0
params['a'].max = 1e6


# Run the fit with Model.fit(Data_Array, Parameters, independent vars).
result = model.fit(S_L1, params, t=t, p_p=p_p)

Upvotes: 1

Views: 1314

Answers (2)

M Newville
M Newville

Reputation: 7862

Your question is stated "how do i write the final fit model (result) expression?". You have answered this question yourself with

def mod_m(t, p_p, a=1, alpha=1,beta=1): 
    return a*(t**alpha)*(p_p**beta)

model = Model(mod_m, independent_vars=['t','p_p'] )

Yes, that is exactly how to write the fitting model.

That, by itself, does not cause the exception

ValueError: The input contains nan values

What causes that ValueError is that your fit function generates nan values with the values of the parameters and independent variables you give it. So... what values are you passing for those?

I recommend printing out values for the parameters in your model function, and values for your independent variables. To be clear, it is easy for exponentiation to produce values larger than 1e308, which will give an inf, and cause the exception you see. So you may have to be more careful about what parameter values are permitted, and this may be sensitive to the values of the independent variables.

Upvotes: 0

James Phillips
James Phillips

Reputation: 4647

Here is a Python 3 example using your function with test data. This uses scipy.optimize.curve_fit() for the multiple regression and creates a 3D data scatterplot, a 3D surface plot of the fitted function, and contour plot of the fitted function. Note that I use the default scipy initial parameters for curve_fit.

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)
    axes = f.add_subplot(111)

    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):
    t = data[0]
    p_p = data[1]
    return a * (t**alpha) * (p_p**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]

    # this example uses curve_fit()'s default initial paramter values
    fittedParameters, pcov = scipy.optimize.curve_fit(func, [xData, yData], zData)

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

    print('fitted prameters', fittedParameters)

Upvotes: 1

Related Questions