Solved – How to run a regression analysis in python on a nonlinear, strange dataset

curve fittingextrapolationpythonregression

I can't figure out what type of regression analysis or extrapolation technique to use in order to come up with an equation for the data I have plotted.

For a school project, I've been testing diodes by connecting them to a temperature switch that shuts off the circuit once the temperature reaches around 250 degrees Celsius. I plotted the data, and I wanted to extrapolate it and get an equation that I could use to predict future diode temperatures. However, I have tried all sorts of extrapolation techniques and none of them accurately fits the data. For example

extrapolation1

as well as this:

extrapolation2

I've tried extrapolation/regressions in both R Studio and python, but I simply cannot get an accurate enough equation for either. Is there anything I could improve in my code, or maybe is there another technique I could use?

This is the code I used to get the first graph:

import numpy as np
import pandas as pd
from scipy.optimize import curve_fit

df = pd.read_csv("JB WELD DIODE REG.csv")
def func(x, a, b, c,d):
    return a*np.sin(b*x-c)+d
x = (df.time)/10**2
y= (df.diode)/10**2


popt, pcov = curve_fit(func, x, y)
popt


plt.plot(x, func(x, *popt), 'r-',
        label='fit: a=%5.3f, b=%5.3f, c=%5.3f, d=%5.3f' %tuple(popt)
        )
plt.plot(x, y)

Best Answer

There are clearly two different segments to each curve. I suggest fitting each separately.

While it is heating up, it looks exponential - see whether log(250-t) is linear.

While it is cooling down, it looks quadratic, so I suggest seeing whether sqrt(t-75) is linear.

Related Solutions

Solved – Strange outcome when performing nonlinear least squares fit to a power law

It would be fairly unusual to find power law data for which the variance is actually constant across $x$ values. However, let's take that constant variance as given and proceed with the analysis as a least squares problem.

The easy way to fit this particular model is actually to use GLMs, not NLS! No need for starting values.

I don't have Matlab on this machine, but in R:

x <- seq(4,6.5,0.5)
y <- 1.0e-04 * c(0.011929047861902, 0.026053683604530, 0.057223759612162,
    0.117413370572612, 0.242357468772138, 0.462327165034928)

 powfit <- glm(y~log(x),family=gaussian(link=log))
 summary(powfit)

[snip]

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -25.05254    0.16275 -153.93 1.07e-08 ***
log(x)        8.05095    0.08824   91.24 8.65e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 6.861052e-14)

    Null deviance: 1.5012e-09  on 5  degrees of freedom
Residual deviance: 2.2273e-13  on 4  degrees of freedom
AIC: -162.52

Number of Fisher Scoring iterations: 1

> plot(x,y,log="xy")
> lines(x,fitted(powfit))

power fit

I highly recommend using Gaussian GLMs with the relevant link function whenever your NLS model can be cast in that form (it can for exponential and power functions, for example).

As to why you're having trouble with it in Matlab, I'm not sure - could it be to do with your convergence criteria? There's no issue with the data at all, but I suggest you start it off at the linear regression coefficients:

> lm(log(y)~log(x))

Call:
lm(formula = log(y) ~ log(x))

Coefficients:
(Intercept)       log(x)  
    -24.203        7.568

So start your parameters at $\exp(-24.2)$ and $7.57$ and see how that works. Or even better, do it the easy way and you can forget about starting values altogether.

Here's my results in R using nls, first, starting with $b$ at the value from the linear log-log fit:

nls(y~a*x^b,start=list(a=exp(-24.2),b=7.57))
Nonlinear regression model
  model: y ~ a * x^b
   data: parent.frame()
        a         b 
1.290e-11 8.063e+00 
 residual sum-of-squares: 2.215e-13

Number of iterations to convergence: 10 
Achieved convergence tolerance: 1.237e-07

And then shifting to starting at b=10:

nls(y~a*x^b,start=list(a=exp(-24.2),b=10))
Nonlinear regression model
  model: y ~ a * x^b
   data: parent.frame()
        a         b 
1.290e-11 8.063e+00 
 residual sum-of-squares: 2.215e-13

Number of iterations to convergence: 19 
Achieved convergence tolerance: 9.267e-08

As you see, they're very close.

The GLM didn't quite match the estimates of the parameters - though a stricter convergence criterion should improve that (glm is only taking one step before concluding it's converged above). It doesn't make much difference to the fit, though. With a bit of fiddling, it's possible to find a closer fit with the GLM (better than the nls), but you start getting down to very tiny deviance (around $10^{-14}$), and I worry if we're really just shifting around numerical error down that far.

Solved – Principal Component Analysis and Regression in Python

Scikit-learn does not have a combined implementation of PCA and regression like for example the pls package in R. But I think one can do like below or choose PLS regression.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

from sklearn.preprocessing import scale
from sklearn.decomposition import PCA
from sklearn import cross_validation
from sklearn.linear_model import LinearRegression

%matplotlib inline

import seaborn as sns
sns.set_style('darkgrid')

df = pd.read_csv('multicollinearity.csv')
X = df.iloc[:,1:6]
y = df.response

Scikit-learn PCA

pca = PCA()

Scale and transform data to get Principal Components

X_reduced = pca.fit_transform(scale(X))

Variance (% cumulative) explained by the principal components

np.cumsum(np.round(pca.explained_variance_ratio_, decimals=4)*100)

array([  73.39,   93.1 ,   98.63,   99.89,  100.  ])

Seems like the first two components indeed explain most of the variance in the data.

10-fold CV, with shuffle

n = len(X_reduced)
kf_10 = cross_validation.KFold(n, n_folds=10, shuffle=True, random_state=2)

regr = LinearRegression()
mse = []

Do one CV to get MSE for just the intercept (no principal components in regression)

score = -1*cross_validation.cross_val_score(regr, np.ones((n,1)), y.ravel(), cv=kf_10, scoring='mean_squared_error').mean()    
mse.append(score)

Do CV for the 5 principle components, adding one component to the regression at the time

for i in np.arange(1,6):
    score = -1*cross_validation.cross_val_score(regr, X_reduced[:,:i], y.ravel(), cv=kf_10, scoring='mean_squared_error').mean()
    mse.append(score)

fig, (ax1, ax2) = plt.subplots(1,2, figsize=(12,5))
ax1.plot(mse, '-v')
ax2.plot([1,2,3,4,5], mse[1:6], '-v')
ax2.set_title('Intercept excluded from plot')

for ax in fig.axes:
    ax.set_xlabel('Number of principal components in regression')
    ax.set_ylabel('MSE')
    ax.set_xlim((-0.2,5.2))

Scikit-learn PLS regression

mse = []

kf_10 = cross_validation.KFold(n, n_folds=10, shuffle=True, random_state=2)

for i in np.arange(1, 6):
    pls = PLSRegression(n_components=i, scale=False)
    pls.fit(scale(X_reduced),y)
    score = cross_validation.cross_val_score(pls, X_reduced, y, cv=kf_10, scoring='mean_squared_error').mean()
    mse.append(-score)

plt.plot(np.arange(1, 6), np.array(mse), '-v')
plt.xlabel('Number of principal components in PLS regression')
plt.ylabel('MSE')
plt.xlim((-0.2, 5.2))