### What is regression?

The purpose of regression is to predict the numerical target value. The most direct way is to write out a calculation formula of the target value according to the input. If you want to predict the power of your little sister's boyfriend's car, you may calculate as follows:

HorsePower = 0.0015 * annualSalary - 0.99 * hoursListeningToPublicRadio

Written in Chinese is:

Car power of little sister's boyfriend = 0.0015 * annual salary of little sister's boyfriend - 0.99 * time to listen to public radio

This is the so-called regression equation, in which 0.0015 and -0.99 are called regression weights. The process of finding these regression coefficients is regression. Once we have these regression coefficients, given the input, it is very easy to make prediction. The specific method is to multiply the regression coefficient by the input value, and then add all the results together to get the predicted value.

When it comes to regression, it generally refers to linear regression, so regression and linear regression in this paper represent the same meaning. Linear regression means that you can multiply the input items by some constants, and then add up the results to get the output. It should be noted that there is another regression model called nonlinear regression, which does not agree with the above practice. For example, it thinks that the output may be the product of the input. In this way, the above power calculation formula can also be written as follows:

HorsePower = 0.0015 * annualSalary / hoursListeningToPublicRadio

This is an example of nonlinear regression, which is not discussed in depth in this paper.

### 1, Find the best fitting line by linear regression

How do we get the regression equation from a lot of data? Suppose that the input data is stored in matrix X and the result is stored in vector y:

The regression coefficients are stored in the vector w:

Then, for the given data x1, that is, the first column data of matrix X, the prediction result u1 will be given by the following formula:

Now the question is, with the data matrix X and the corresponding label vector Y in hand, how can we find w? A common method is to find the w that minimizes the error. The error here refers to the difference between the predicted u value and the real y value. The simple accumulation of the error will make the positive difference and negative difference offset each other, so we use the square error.

The sum of squared errors can be written as:

The matrix representation can also be written as:

Why can it change like this? Remember a premise: if x is a vector, X is the column vector and x^T is the row vector by default. Bring in the data matrix X and label vector y mentioned above to know why it changes so much.

Before we continue the derivation, we should first clarify one purpose: find w and minimize the sum of square errors. Because we believe that the smaller the sum of square errors, the better the fitting effect of linear regression.

Now, we take the derivative of w by the sum of the square errors expressed by the matrix:

If you are not familiar with matrix calculation, you can move here: https://blog.csdn.net/nomadlx53/article/details/50849941

Let the above formula be equal to 0 to obtain:

The small mark above w indicates that this is the optimal solution of w that can be estimated at present. The w estimated from the existing data may not be the real W value in the data, so a "hat" symbol is used here to indicate that it is only a best estimate of W.

It is worth noting that the above formula contains an inverse matrix, that is, the equation is only used when the inverse matrix exists, that is, the matrix is a square matrix and its determinant is not 0.

The best w solution is a common problem in statistics. There are many other methods besides matrix method. By calling the matrix method in NumPy library, we can complete the required function with only a few lines of code. This method is also called OLS, which means "ordinary least squares".

Let's look at the dataset first:

The first column is 1.0, i.e. x0. The second column is x1, i.e. X-axis data. The third column is x2, i.e. y-axis data. First, draw the data and look at the data distribution. The code is as follows:

import matplotlib.pyplot as plt import numpy as np def loadDataSet(fileName): """ Function description:Load data Parameters: fileName - file name Returns: xArr - x data set yArr - y data set """ numFeat = len(open(fileName).readline().split('\t')) - 1 xArr = []; yArr = [] fr = open(fileName) for line in fr.readlines(): lineArr =[] curLine = line.strip().split('\t') for i in range(numFeat): lineArr.append(float(curLine[i])) xArr.append(lineArr) yArr.append(float(curLine[-1])) return xArr, yArr def plotDataSet(): """ Function description:Draw dataset Parameters: nothing Returns: nothing """ xArr, yArr = loadDataSet('ex0.txt') #Load dataset n = len(xArr) #Number of data xcord = []; ycord = [] #Sample point for i in range(n): xcord.append(xArr[i][1]); ycord.append(yArr[i]) #Sample point fig = plt.figure() ax = fig.add_subplot(111) #Add subplot ax.scatter(xcord, ycord, s = 20, c = 'blue',alpha = .5) #Draw sample points plt.title('DataSet') #Draw title plt.xlabel('X') plt.show() if __name__ == '__main__': plotDataSet()

By visualizing the data, we can see the distribution of the data. Next, let's calculate the regression coefficient vector according to the regression coefficient calculation method derived above, and draw the regression curve according to the regression coefficient vector. The code is as follows:

import matplotlib.pyplot as plt import numpy as np def loadDataSet(fileName): """ Function description:Load data Parameters: fileName - file name Returns: xArr - x data set yArr - y data set """ numFeat = len(open(fileName).readline().split('\t')) - 1 xArr = []; yArr = [] fr = open(fileName) for line in fr.readlines(): lineArr =[] curLine = line.strip().split('\t') for i in range(numFeat): lineArr.append(float(curLine[i])) xArr.append(lineArr) yArr.append(float(curLine[-1])) return xArr, yArr def standRegres(xArr,yArr): """ Function description:Calculate regression coefficient w Parameters: xArr - x data set yArr - y data set Returns: ws - regression coefficient """ xMat = np.mat(xArr); yMat = np.mat(yArr).T xTx = xMat.T * xMat #The regression coefficient is calculated according to the publicity derived in this paper if np.linalg.det(xTx) == 0.0: print("The matrix is singular,No inversion") return ws = xTx.I * (xMat.T*yMat) return ws def plotRegression(): """ Function description:Draw regression curves and data points Parameters: nothing Returns: nothing """ xArr, yArr = loadDataSet('ex0.txt') #Load dataset ws = standRegres(xArr, yArr) #Calculate regression coefficient xMat = np.mat(xArr) #Create xMat matrix yMat = np.mat(yArr) #Create yMat matrix xCopy = xMat.copy() #Deep copy xMat matrix xCopy.sort(0) #sort yHat = xCopy * ws #Calculate the corresponding y value fig = plt.figure() ax = fig.add_subplot(111) #Add subplot ax.plot(xCopy[:, 1], yHat, c = 'red') #Draw regression curve ax.scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = 'blue',alpha = .5) #Draw sample points plt.title('DataSet') #Draw title plt.xlabel('X') plt.show() if __name__ == '__main__': plotRegression()

How to judge the fitting effect of the fitting curve? Of course, we can observe according to our own experience. In addition, we can also use the corrcoef method to compare the correlation between the predicted value and the real value. The code is as follows:

import numpy as np def loadDataSet(fileName): """ Function description:Load data Parameters: fileName - file name Returns: xArr - x data set yArr - y data set """ numFeat = len(open(fileName).readline().split('\t')) - 1 xArr = []; yArr = [] fr = open(fileName) for line in fr.readlines(): lineArr =[] curLine = line.strip().split('\t') for i in range(numFeat): lineArr.append(float(curLine[i])) xArr.append(lineArr) yArr.append(float(curLine[-1])) return xArr, yArr def standRegres(xArr,yArr): """ Function description:Calculate regression coefficient w Parameters: xArr - x data set yArr - y data set Returns: ws - regression coefficient """ xMat = np.mat(xArr); yMat = np.mat(yArr).T xTx = xMat.T * xMat #The regression coefficient is calculated according to the publicity derived in this paper if np.linalg.det(xTx) == 0.0: print("The matrix is singular,No inversion") return ws = xTx.I * (xMat.T*yMat) return ws if __name__ == '__main__': xArr, yArr = loadDataSet('ex0.txt') #Load dataset ws = standRegres(xArr, yArr) #Calculate regression coefficient xMat = np.mat(xArr) #Create xMat matrix yMat = np.mat(yArr) #Create yMat matrix yHat = xMat * ws print(np.corrcoef(yHat.T, yMat))

It can be seen that the data on the diagonal is 1.0, because the match between yMat and itself is perfect, while the correlation coefficient between yhhat and yMat is 0.98.

The best fitting method regards the data as a straight line for modeling, which has a very good performance. There seem to be other potential patterns in the data. So how can we use these patterns? We can locally adjust the forecast according to the data. This method will be introduced below.

### 2, Locally weighted linear regression

One problem of linear regression is that under fitting may occur, because it seeks unbiased estimation with small mean square error. Obviously, if the model is under fitted, it will not achieve good prediction results. Therefore, some methods allow some deviations to be introduced into the estimation, so as to reduce the mean square error of prediction.

One method is Locally Weighted Linear Regression (LWLR). In this method, we give a certain weight to each point near the point to be predicted. Like kNN, this algorithm needs to select the corresponding data subset in advance for each prediction. The form of the regression coefficient W released by the algorithm is as follows:

Where W is a matrix. The difference between this formula and the formula derived above is that w is used to give weight to each store.

LWLR uses a "kernel" (similar to the kernel in support vector machine) to give higher weights to nearby points. The type of kernel can be selected freely. The most commonly used kernel is Gaussian kernel. The weight of Gaussian kernel is as follows:

In this way, we can write locally weighted linear regression according to the above formula. We can adjust the regression effect by changing the value of k. The code is as follows:

from matplotlib.font_manager import FontProperties import matplotlib.pyplot as plt import numpy as np def loadDataSet(fileName): """ Function description:Load data Parameters: fileName - file name Returns: xArr - x data set yArr - y data set """ numFeat = len(open(fileName).readline().split('\t')) - 1 xArr = []; yArr = [] fr = open(fileName) for line in fr.readlines(): lineArr =[] curLine = line.strip().split('\t') for i in range(numFeat): lineArr.append(float(curLine[i])) xArr.append(lineArr) yArr.append(float(curLine[-1])) return xArr, yArr def plotlwlrRegression(): """ Function description:Draw multiple locally weighted regression curves Parameters: nothing Returns: nothing """ font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14) xArr, yArr = loadDataSet('ex0.txt') #Load dataset yHat_1 = lwlrTest(xArr, xArr, yArr, 1.0) #Yhhat was calculated from locally weighted linear regression yHat_2 = lwlrTest(xArr, xArr, yArr, 0.01) #Yhhat was calculated from locally weighted linear regression yHat_3 = lwlrTest(xArr, xArr, yArr, 0.003) #Yhhat was calculated from locally weighted linear regression xMat = np.mat(xArr) #Create xMat matrix yMat = np.mat(yArr) #Create yMat matrix srtInd = xMat[:, 1].argsort(0) #Sort, return index value xSort = xMat[srtInd][:,0,:] fig, axs = plt.subplots(nrows=3, ncols=1,sharex=False, sharey=False, figsize=(10,8)) axs[0].plot(xSort[:, 1], yHat_1[srtInd], c = 'red') #Draw regression curve axs[1].plot(xSort[:, 1], yHat_2[srtInd], c = 'red') #Draw regression curve axs[2].plot(xSort[:, 1], yHat_3[srtInd], c = 'red') #Draw regression curve axs[0].scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = 'blue', alpha = .5) #Draw sample points axs[1].scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = 'blue', alpha = .5) #Draw sample points axs[2].scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = 'blue', alpha = .5) #Draw sample points #Set title, x-axis label,y-axis label axs0_title_text = axs[0].set_title(u'Locally weighted regression curve,k=1.0',FontProperties=font) axs1_title_text = axs[1].set_title(u'Locally weighted regression curve,k=0.01',FontProperties=font) axs2_title_text = axs[2].set_title(u'Locally weighted regression curve,k=0.003',FontProperties=font) plt.setp(axs0_title_text, size=8, weight='bold', color='red') plt.setp(axs1_title_text, size=8, weight='bold', color='red') plt.setp(axs2_title_text, size=8, weight='bold', color='red') plt.xlabel('X') plt.show() def lwlr(testPoint, xArr, yArr, k = 1.0): """ Function description:Regression coefficients were calculated using locally weighted linear regression w Parameters: testPoint - Test sample point xArr - x data set yArr - y data set k - Gaussian kernel k,Custom parameters Returns: ws - regression coefficient """ xMat = np.mat(xArr); yMat = np.mat(yArr).T m = np.shape(xMat)[0] weights = np.mat(np.eye((m))) #Create weight diagonal matrix for j in range(m): #Traverse the data set to calculate the weight of each sample diffMat = testPoint - xMat[j, :] weights[j, j] = np.exp(diffMat * diffMat.T/(-2.0 * k**2)) xTx = xMat.T * (weights * xMat) if np.linalg.det(xTx) == 0.0: print("The matrix is singular,No inversion") return ws = xTx.I * (xMat.T * (weights * yMat)) #Calculate regression coefficient return testPoint * ws def lwlrTest(testArr, xArr, yArr, k=1.0): """ Function description:Locally weighted linear regression test Parameters: testArr - Test data set xArr - x data set yArr - y data set k - Gaussian kernel k,Custom parameters Returns: ws - regression coefficient """ m = np.shape(testArr)[0] #Calculate test dataset size yHat = np.zeros(m) for i in range(m): #Predict each sample point yHat[i] = lwlr(testArr[i],xArr,yArr,k) return yHat if __name__ == '__main__': plotlwlrRegression()

It can be seen that the smaller the K, the better the fitting effect. However, when k is too small, over fitting occurs, for example, when k is equal to 0.003.

### 3, Example: predicting the age of abalone

Next, we apply regression to real data. abalone.txt file records the age of Abalone (an aquatic organism → _→), which is from the data of UCI data set. Abalone age can be calculated from the number of layers of abalone shell.

The data of the dataset is as follows:

As you can see, the data set is multidimensional, so it is difficult for us to draw its distribution. The meaning of the representation of each dimension data is not given, but it doesn't matter. We just need to know that the data in the last column is the y value. The last column represents the real age of abalone. The data in the previous columns are the characteristics of abalone, such as the number of layers of abalone shell. Without data cleaning, we directly use all features to test our locally weighted regression.

Create a new. py file and add rssError function to evaluate the final regression results. The code is as follows:

from matplotlib.font_manager import FontProperties import matplotlib.pyplot as plt import numpy as np def loadDataSet(fileName): """ Function description:Load data Parameters: fileName - file name Returns: xArr - x data set yArr - y data set """ numFeat = len(open(fileName).readline().split('\t')) - 1 xArr = []; yArr = [] fr = open(fileName) for line in fr.readlines(): lineArr =[] curLine = line.strip().split('\t') for i in range(numFeat): lineArr.append(float(curLine[i])) xArr.append(lineArr) yArr.append(float(curLine[-1])) return xArr, yArr def lwlr(testPoint, xArr, yArr, k = 1.0): """ Function description:Regression coefficients were calculated using locally weighted linear regression w Parameters: testPoint - Test sample point xArr - x data set yArr - y data set k - Gaussian kernel k,Custom parameters Returns: ws - regression coefficient """ xMat = np.mat(xArr); yMat = np.mat(yArr).T m = np.shape(xMat)[0] weights = np.mat(np.eye((m))) #Create weight diagonal matrix for j in range(m): #Traverse the data set to calculate the weight of each sample diffMat = testPoint - xMat[j, :] weights[j, j] = np.exp(diffMat * diffMat.T/(-2.0 * k**2)) xTx = xMat.T * (weights * xMat) if np.linalg.det(xTx) == 0.0: print("The matrix is singular,No inversion") return ws = xTx.I * (xMat.T * (weights * yMat)) #Calculate regression coefficient return testPoint * ws def lwlrTest(testArr, xArr, yArr, k=1.0): """ Function description:Locally weighted linear regression test Parameters: testArr - Test data set,Test set xArr - x data set,Training set yArr - y data set,Training set k - Gaussian kernel k,Custom parameters Returns: ws - regression coefficient """ m = np.shape(testArr)[0] #Calculate test dataset size yHat = np.zeros(m) for i in range(m): #Predict each sample point yHat[i] = lwlr(testArr[i],xArr,yArr,k) return yHat def standRegres(xArr,yArr): """ Function description:Calculate regression coefficient w Parameters: xArr - x data set yArr - y data set Returns: ws - regression coefficient """ xMat = np.mat(xArr); yMat = np.mat(yArr).T xTx = xMat.T * xMat #The regression coefficient is calculated according to the publicity derived in this paper if np.linalg.det(xTx) == 0.0: print("The matrix is singular,No inversion") return ws = xTx.I * (xMat.T*yMat) return ws def rssError(yArr, yHatArr): """ Error size evaluation function Parameters: yArr - Real data yHatArr - Forecast data Returns: Error size """ return ((yArr - yHatArr) **2).sum() if __name__ == '__main__': abX, abY = loadDataSet('abalone.txt') print('The training set is the same as the test set:Locally weighted linear regression,nucleus k Effect of the size of on the prediction:') yHat01 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 0.1) yHat1 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 1) yHat10 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 10) print('k=0.1 Time,The error size is:',rssError(abY[0:99], yHat01.T)) print('k=1 Time,The error size is:',rssError(abY[0:99], yHat1.T)) print('k=10 Time,The error size is:',rssError(abY[0:99], yHat10.T)) print('') print('The training set is different from the test set:Locally weighted linear regression,nucleus k Is the smaller the better? Replace dataset,The test results are as follows:') yHat01 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 0.1) yHat1 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 1) yHat10 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 10) print('k=0.1 Time,The error size is:',rssError(abY[100:199], yHat01.T)) print('k=1 Time,The error size is:',rssError(abY[100:199], yHat1.T)) print('k=10 Time,The error size is:',rssError(abY[100:199], yHat10.T)) print('') print('The training set is different from the test set:Simple linear regression and k=1 Comparison of locally weighted linear regression:') print('k=1 Time,The error size is:', rssError(abY[100:199], yHat1.T)) ws = standRegres(abX[0:99], abY[0:99]) yHat = np.mat(abX[100:199]) * ws print('Simple linear regression error size:', rssError(abY[100:199], yHat.T.A))

It can be seen that when k=0.1, the error of the training set is small, but when applied to the new data set, the error becomes larger. This is often called over fitting phenomenon. For the model we train, we should ensure the high accuracy of the test set, so that the trained model can be applied to new data, that is, we should strengthen the universality of the model. It can be seen that when k=1, the effect of locally weighted linear regression is similar to that of simple linear regression. This also shows that the best model can only be selected by comparing the effects on unknown data. So the best nuclear size is 10? Maybe, but if you want to get better results, you should do 10 tests with 10 different sample sets to compare the results.

This example shows how to use locally weighted linear regression to build a model, which can get better results than ordinary linear regression. The problem with locally weighted linear regression is that it must be run on the whole data set each time. In other words, in order to make a prediction, all training data must be saved.

### 4, Reduction factor to "understand" data

#### 1. Ridge regression

from matplotlib.font_manager import FontProperties import matplotlib.pyplot as plt import numpy as np def loadDataSet(fileName): """ Function description:Load data Parameters: fileName - file name Returns: xArr - x data set yArr - y data set """ numFeat = len(open(fileName).readline().split('\t')) - 1 xArr = []; yArr = [] fr = open(fileName) for line in fr.readlines(): lineArr =[] curLine = line.strip().split('\t') for i in range(numFeat): lineArr.append(float(curLine[i])) xArr.append(lineArr) yArr.append(float(curLine[-1])) return xArr, yArr def ridgeRegres(xMat, yMat, lam = 0.2): """ Function description:Ridge regression Parameters: xMat - x data set yMat - y data set lam - Reduction factor Returns: ws - regression coefficient """ xTx = xMat.T * xMat denom = xTx + np.eye(np.shape(xMat)[1]) * lam if np.linalg.det(denom) == 0.0: print("The matrix is singular,Cannot transpose") return ws = denom.I * (xMat.T * yMat) return ws def ridgeTest(xArr, yArr): """ Function description:Ridge regression test Parameters: xMat - x data set yMat - y data set Returns: wMat - Regression coefficient matrix """ xMat = np.mat(xArr); yMat = np.mat(yArr).T #Data standardization yMean = np.mean(yMat, axis = 0) #Row to row operation, find the mean yMat = yMat - yMean #Data minus mean xMeans = np.mean(xMat, axis = 0) #Row to row operation, find the mean xVar = np.var(xMat, axis = 0) #Row to row operation, variance xMat = (xMat - xMeans) / xVar #Data is normalized by subtracting mean divided by variance numTestPts = 30 #30 different lambda tests wMat = np.zeros((numTestPts, np.shape(xMat)[1])) #Initial regression coefficient matrix for i in range(numTestPts): #Changing lambda to calculate regression coefficient ws = ridgeRegres(xMat, yMat, np.exp(i - 10)) #lambda changes exponentially with e, initially a very small number, wMat[i, :] = ws.T #Calculate regression coefficient matrix return wMat def plotwMat(): """ Function description:Plot ridge regression coefficient matrix """ font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14) abX, abY = loadDataSet('abalone.txt') redgeWeights = ridgeTest(abX, abY) fig = plt.figure() ax = fig.add_subplot(111) ax.plot(redgeWeights) ax_title_text = ax.set_title(u'log(lambada)Relationship with regression coefficient', FontProperties = font) ax_xlabel_text = ax.set_xlabel(u'log(lambada)', FontProperties = font) ax_ylabel_text = ax.set_ylabel(u'regression coefficient', FontProperties = font) plt.setp(ax_title_text, size = 20, weight = 'bold', color = 'red') plt.setp(ax_xlabel_text, size = 10, weight = 'bold', color = 'black') plt.setp(ax_ylabel_text, size = 10, weight = 'bold', color = 'black') plt.show() if __name__ == '__main__': plotwMat()

The regression coefficients and log are plotted( λ) On the far left, i.e λ At the minimum, the original values of all coefficients can be obtained (consistent with linear regression); on the right, the coefficients are reduced to 0; at a certain position in the middle, the best prediction result will be obtained. Want to get the best prediction result λ Parameters can be obtained by cross validation, which will be explained later in the article.

#### 2,lasso

It is not difficult to prove that when the following constraints are added, the ordinary least squares regression will obtain the same formula as the ridge regression:

The above formula defines that the sum of squares of all regression coefficients cannot be greater than λ. When two or more features are correlated, a large positive coefficient and a large negative number may be obtained by ordinary least squares regression. Because of the above constraints, ridge regression can avoid this problem.

Similar to ridge regression, another reduction method lasso also limits the regression coefficient, and the corresponding constraints are as follows:

The only difference is that this constraint replaces the sum of squares with absolute values. Although the form of constraints changes only slightly, the results are quite different: in λ When it is small enough, some coefficients will be forced to reduce to 0. This feature can help us better understand the data. The two constraints look similar in terms of formula, but slight changes greatly increase the computational complexity (in order to solve the regression coefficient under this new constraint, the quadratic programming algorithm needs to be used). The following will introduce a simpler method to obtain the results, which is called forward stepwise regression.

#### 3. Forward stepwise regression

The forward stepwise linear regression algorithm belongs to a greedy algorithm, that is, the error is reduced as much as possible at each step. We calculate the regression coefficient no longer by formula, but by fine tuning each regression coefficient every time, and then calculating the prediction error. The set of regression coefficients that minimize the error is the best regression coefficient we need.

The implementation of forward stepwise linear regression is also very simple. Of course, the data should be standardized first. The code is as follows:

from matplotlib.font_manager import FontProperties import matplotlib.pyplot as plt import numpy as np def loadDataSet(fileName): """ Function description:Load data Parameters: fileName - file name Returns: xArr - x data set yArr - y data set """ numFeat = len(open(fileName).readline().split('\t')) - 1 xArr = []; yArr = [] fr = open(fileName) for line in fr.readlines(): lineArr =[] curLine = line.strip().split('\t') for i in range(numFeat): lineArr.append(float(curLine[i])) xArr.append(lineArr) yArr.append(float(curLine[-1])) return xArr, yArr def regularize(xMat, yMat): """ Function description:Data standardization Parameters: xMat - x data set yMat - y data set Returns: inxMat - Standardized x data set inyMat - Standardized y data set """ inxMat = xMat.copy() #Data copy inyMat = yMat.copy() yMean = np.mean(yMat, 0) #Row to row operation, find the mean inyMat = yMat - yMean #Data minus mean inMeans = np.mean(inxMat, 0) #Row to row operation, find the mean inVar = np.var(inxMat, 0) #Row to row operation, variance inxMat = (inxMat - inMeans) / inVar #Data is normalized by subtracting mean divided by variance return inxMat, inyMat def rssError(yArr,yHatArr): """ Function description:Calculate square error Parameters: yArr - Estimate yHatArr - True value Returns: """ return ((yArr-yHatArr)**2).sum() def stageWise(xArr, yArr, eps = 0.01, numIt = 100): """ Function description:Forward stepwise linear regression Parameters: xArr - x input data yArr - y Forecast data eps - Step size to be adjusted for each iteration numIt - Number of iterations Returns: returnMat - numIt Regression coefficient matrix of second iteration """ xMat = np.mat(xArr); yMat = np.mat(yArr).T #data set xMat, yMat = regularize(xMat, yMat) #Data standardization m, n = np.shape(xMat) returnMat = np.zeros((numIt, n)) #Initialize the regression coefficient matrix of numIt iterations ws = np.zeros((n, 1)) #Initialize regression coefficient matrix wsTest = ws.copy() wsMax = ws.copy() for i in range(numIt): #Iteration numIt times # print(ws.T) #Print current regression coefficient matrix lowestError = float('inf'); #Positive infinity for j in range(n): #Traverse the regression coefficients of each feature for sign in [-1, 1]: wsTest = ws.copy() wsTest[j] += eps * sign #Fine tuning regression coefficient yTest = xMat * wsTest #Calculate predicted value rssE = rssError(yMat.A, yTest.A) #Calculate square error if rssE < lowestError: #If the error is smaller, the current optimal regression coefficient is updated lowestError = rssE wsMax = wsTest ws = wsMax.copy() returnMat[i,:] = ws.T #Record the regression coefficient matrix of numIt iterations return returnMat def plotstageWiseMat(): """ Function description:Plot ridge regression coefficient matrix """ font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14) xArr, yArr = loadDataSet('abalone.txt') returnMat = stageWise(xArr, yArr, 0.005, 1000) fig = plt.figure() ax = fig.add_subplot(111) ax.plot(returnMat) ax_title_text = ax.set_title(u'Forward stepwise regression:Relationship between iteration number and regression coefficient', FontProperties = font) ax_xlabel_text = ax.set_xlabel(u'Number of iterations', FontProperties = font) ax_ylabel_text = ax.set_ylabel(u'regression coefficient', FontProperties = font) plt.setp(ax_title_text, size = 15, weight = 'bold', color = 'red') plt.setp(ax_xlabel_text, size = 10, weight = 'bold', color = 'black') plt.setp(ax_ylabel_text, size = 10, weight = 'bold', color = 'black') plt.show() if __name__ == '__main__': plotstageWiseMat()

Still, we printed the relationship curve between the number of iterations and the regression coefficient. It can be seen that some coefficients are about 0 from beginning to end, which shows that they do not have any impact on the target, that is, these characteristics are probably unnecessary. The advantage of the stepwise linear regression algorithm is that it can help people understand some models and make improvements. After building a model , the algorithm can be run to find out the important features, so it is possible to stop the collection of those unimportant features in time.

To sum up:

The reduction method (stepwise linear regression or ridge regression) is to reduce some coefficients to very small values or directly to 0. This will increase the deviation of the model (reduce the weight of some features) , by reducing the regression coefficient of some features to 0, the complexity of the model is reduced. After eliminating the redundant features, the model is easier to understand and the prediction error is reduced. However, when the reduction is too strict, there will be over fitting, that is, the prediction result with the training set is very good, and the prediction with the test set is much worse.

### 5, Example: forecasting the price of LEGO toy sets

LEGO produces assembled toys, which are composed of many plastic inserts of different sizes. Its appearance is shown in the figure below:

Generally speaking, these inserts are sold in sets. They can be assembled into many different things, such as ships, castles, some famous buildings, etc. the number of parts contained in each set of LEGO company ranges from 10 to 5000.

A Lego kit will basically be discontinued in a few years, but LEGO collectors will still trade with each other after shutdown. This example is to use the regression method to predict the transaction price between collectors.

#### 1. Get data

We get the information we need by parsing the html file:

from bs4 import BeautifulSoup def scrapePage(retX, retY, inFile, yr, numPce, origPrc): """ Function description:Read data from the page and generate retX and retY list Parameters: retX - data X retY - data Y inFile - HTML file yr - particular year numPce - Number of LEGO parts origPrc - original price Returns: nothing """ # Open and read HTML file with open(inFile, encoding='utf-8') as f: html = f.read() soup = BeautifulSoup(html) i = 1 # Parse according to HTML page structure currentRow = soup.find_all('table', r = "%d" % i) while(len(currentRow) != 0): currentRow = soup.find_all('table', r = "%d" % i) title = currentRow[0].find_all('a')[1].text lwrTitle = title.lower() # Find out if there are brand new labels if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1): newFlag = 1.0 else: newFlag = 0.0 # To find out whether the logo has been sold, we only collect the data that has been sold soldUnicde = currentRow[0].find_all('td')[3].find_all('span') if len(soldUnicde) == 0: print("commodity #%d not sold "% i) else: # Analyze the page to get the current price soldPrice = currentRow[0].find_all('td')[4] priceStr = soldPrice.text priceStr = priceStr.replace('$','') priceStr = priceStr.replace(',','') if len(soldPrice) > 1: priceStr = priceStr.replace('Free shipping', '') sellingPrice = float(priceStr) # Remove incomplete package price if sellingPrice > origPrc * 0.5: print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice)) retX.append([yr, numPce, newFlag, origPrc]) retY.append(sellingPrice) i += 1 currentRow = soup.find_all('table', r = "%d" % i) def setDataCollect(retX, retY): """ Function description:Read the data of six Lego sets in turn and generate the data matrix Parameters: nothing Returns: nothing """ scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99) #LEGO 8288 in 2006, number of parts 800, original price 49.99 scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99) #LEGO 10030 in 2002, number of parts 3096, original price 269.99 scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99) #LEGO 10179 in 2007, number of parts 5195, original price 499.99 scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99) #LEGO 10181 in 2007, number of parts 3428, original price 199.99 scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99) #LEGO 10189 in 2008, number of parts 5922, original price 299.99 scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99) #LEGO 10196 in 2009, 3263 parts, original price 249.99 if __name__ == '__main__': lgX = [] lgY = [] setDataCollect(lgX, lgY)

We deal with goods that are not available. These characteristics are: year of production, number of parts, whether it is brand new, original price and selling price (second-hand transaction).

#### 2. Model building

After we have processed the data set, the next step is to train the model. First, we need to add the feature X0 column with all 0. Because the first column of linear regression requires 1.0. Then, use the simplest ordinary linear regression i, and write the code as follows:

import numpy as np from bs4 import BeautifulSoup def scrapePage(retX, retY, inFile, yr, numPce, origPrc): """ Function description:Read data from the page and generate retX and retY list Parameters: retX - data X retY - data Y inFile - HTML file yr - particular year numPce - Number of LEGO parts origPrc - original price Returns: nothing """ # Open and read HTML file with open(inFile, encoding='utf-8') as f: html = f.read() soup = BeautifulSoup(html) i = 1 # Parse according to HTML page structure currentRow = soup.find_all('table', r = "%d" % i) while(len(currentRow) != 0): currentRow = soup.find_all('table', r = "%d" % i) title = currentRow[0].find_all('a')[1].text lwrTitle = title.lower() # Find out if there are brand new labels if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1): newFlag = 1.0 else: newFlag = 0.0 # To find out whether the logo has been sold, we only collect the data that has been sold soldUnicde = currentRow[0].find_all('td')[3].find_all('span') if len(soldUnicde) == 0: print("commodity #%d not sold "% i) else: # Analyze the page to get the current price soldPrice = currentRow[0].find_all('td')[4] priceStr = soldPrice.text priceStr = priceStr.replace('$','') priceStr = priceStr.replace(',','') if len(soldPrice) > 1: priceStr = priceStr.replace('Free shipping', '') sellingPrice = float(priceStr) # Remove incomplete package price if sellingPrice > origPrc * 0.5: print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice)) retX.append([yr, numPce, newFlag, origPrc]) retY.append(sellingPrice) i += 1 currentRow = soup.find_all('table', r = "%d" % i) def setDataCollect(retX, retY): """ Function description:Read the data of six Lego sets in turn and generate the data matrix Parameters: nothing Returns: nothing """ scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99) #LEGO 8288 in 2006, number of parts 800, original price 49.99 scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99) #LEGO 10030 in 2002, number of parts 3096, original price 269.99 scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99) #LEGO 10179 in 2007, number of parts 5195, original price 499.99 scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99) #LEGO 10181 in 2007, number of parts 3428, original price 199.99 scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99) #LEGO 10189 in 2008, number of parts 5922, original price 299.99 scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99) #LEGO 10196 in 2009, 3263 parts, original price 249.99 def regularize(xMat, yMat): """ Function description:Data standardization Parameters: xMat - x data set yMat - y data set Returns: inxMat - Standardized x data set inyMat - Standardized y data set """ inxMat = xMat.copy() #Data copy inyMat = yMat.copy() yMean = np.mean(yMat, 0) #Row to row operation, find the mean inyMat = yMat - yMean #Data minus mean inMeans = np.mean(inxMat, 0) #Row to row operation, find the mean inVar = np.var(inxMat, 0) #Row to row operation, variance # print(inxMat) print(inMeans) # print(inVar) inxMat = (inxMat - inMeans) / inVar #Data is normalized by subtracting mean divided by variance return inxMat, inyMat def rssError(yArr,yHatArr): """ Function description:Calculate square error Parameters: yArr - Estimate yHatArr - True value Returns: """ return ((yArr-yHatArr)**2).sum() def standRegres(xArr,yArr): """ Function description:Calculate regression coefficient w Parameters: xArr - x data set yArr - y data set Returns: ws - regression coefficient """ xMat = np.mat(xArr); yMat = np.mat(yArr).T xTx = xMat.T * xMat #The regression coefficient is calculated according to the publicity derived in this paper if np.linalg.det(xTx) == 0.0: print("The matrix is singular,Cannot transpose") return ws = xTx.I * (xMat.T*yMat) return ws def useStandRegres(): """ Function description:Using simple linear regression Parameters: nothing Returns: nothing """ lgX = [] lgY = [] setDataCollect(lgX, lgY) data_num, features_num = np.shape(lgX) lgX1 = np.mat(np.ones((data_num, features_num + 1))) lgX1[:, 1:5] = np.mat(lgX) ws = standRegres(lgX1, lgY) print('%f%+f*particular year%+f*Number of parts%+f*Is it brand new%+f*original price' % (ws[0],ws[1],ws[2],ws[3],ws[4])) if __name__ == '__main__': useStandRegres()

It can be seen that the formula used in the model is shown in the figure above. Although this model fits the data well, it doesn't seem to make any sense. The more parts in the kit, the lower the price, which is unreasonable.

We use ridge regression to find the method to minimize the error through cross validation λ Corresponding regression coefficient. The code is as follows:

import numpy as np from bs4 import BeautifulSoup import random def scrapePage(retX, retY, inFile, yr, numPce, origPrc): """ Function description:Read data from the page and generate retX and retY list Parameters: retX - data X retY - data Y inFile - HTML file yr - particular year numPce - Number of LEGO parts origPrc - original price Returns: nothing """ # Open and read HTML file with open(inFile, encoding='utf-8') as f: html = f.read() soup = BeautifulSoup(html) i = 1 # Parse according to HTML page structure currentRow = soup.find_all('table', r = "%d" % i) while(len(currentRow) != 0): currentRow = soup.find_all('table', r = "%d" % i) title = currentRow[0].find_all('a')[1].text lwrTitle = title.lower() # Find out if there are brand new labels if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1): newFlag = 1.0 else: newFlag = 0.0 # To find out whether the logo has been sold, we only collect the data that has been sold soldUnicde = currentRow[0].find_all('td')[3].find_all('span') if len(soldUnicde) == 0: print("commodity #%d not sold "% i) else: # Analyze the page to get the current price soldPrice = currentRow[0].find_all('td')[4] priceStr = soldPrice.text priceStr = priceStr.replace('$','') priceStr = priceStr.replace(',','') if len(soldPrice) > 1: priceStr = priceStr.replace('Free shipping', '') sellingPrice = float(priceStr) # Remove incomplete package price if sellingPrice > origPrc * 0.5: print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice)) retX.append([yr, numPce, newFlag, origPrc]) retY.append(sellingPrice) i += 1 currentRow = soup.find_all('table', r = "%d" % i) def ridgeRegres(xMat, yMat, lam = 0.2): """ Function description:Ridge regression Parameters: xMat - x data set yMat - y data set lam - Reduction factor Returns: ws - regression coefficient """ xTx = xMat.T * xMat denom = xTx + np.eye(np.shape(xMat)[1]) * lam if np.linalg.det(denom) == 0.0: print("The matrix is singular,Cannot transpose") return ws = denom.I * (xMat.T * yMat) return ws def setDataCollect(retX, retY): """ Function description:Read the data of six Lego sets in turn and generate the data matrix Parameters: nothing Returns: nothing """ scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99) #LEGO 8288 in 2006, number of parts 800, original price 49.99 scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99) #LEGO 10030 in 2002, number of parts 3096, original price 269.99 scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99) #LEGO 10179 in 2007, number of parts 5195, original price 499.99 scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99) #LEGO 10181 in 2007, number of parts 3428, original price 199.99 scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99) #LEGO 10189 in 2008, number of parts 5922, original price 299.99 scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99) #LEGO 10196 in 2009, 3263 parts, original price 249.99 def regularize(xMat, yMat): """ Function description:Data standardization Parameters: xMat - x data set yMat - y data set Returns: inxMat - Standardized x data set inyMat - Standardized y data set """ inxMat = xMat.copy() #Data copy inyMat = yMat.copy() yMean = np.mean(yMat, 0) #Row to row operation, find the mean inyMat = yMat - yMean #Data minus mean inMeans = np.mean(inxMat, 0) #Row to row operation, find the mean inVar = np.var(inxMat, 0) #Row to row operation, variance # print(inxMat) print(inMeans) # print(inVar) inxMat = (inxMat - inMeans) / inVar #Data is normalized by subtracting mean divided by variance return inxMat, inyMat def rssError(yArr,yHatArr): """ Function description:Calculate square error Parameters: yArr - Estimate yHatArr - True value Returns: """ return ((yArr-yHatArr)**2).sum() def standRegres(xArr,yArr): """ Function description:Calculate regression coefficient w Parameters: xArr - x data set yArr - y data set Returns: ws - regression coefficient """ xMat = np.mat(xArr); yMat = np.mat(yArr).T xTx = xMat.T * xMat #The regression coefficient is calculated according to the publicity derived in this paper if np.linalg.det(xTx) == 0.0: print("The matrix is singular,Cannot transpose") return ws = xTx.I * (xMat.T*yMat) return ws def crossValidation(xArr, yArr, numVal = 10): """ Function description:Cross validation ridge regression Parameters: xArr - x data set yArr - y data set numVal - Number of cross validation Returns: wMat - Regression coefficient matrix """ m = len(yArr) #Number of statistical samples indexList = list(range(m)) #Generate index value list errorMat = np.zeros((numVal,30)) #create error mat 30columns numVal rows for i in range(numVal): #Cross validation numVal times trainX = []; trainY = [] #Training set testX = []; testY = [] #Test set random.shuffle(indexList) #Disorganize for j in range(m): #Divided data set: 90% training set and 10% test set if j < m * 0.9: trainX.append(xArr[indexList[j]]) trainY.append(yArr[indexList[j]]) else: testX.append(xArr[indexList[j]]) testY.append(yArr[indexList[j]]) wMat = ridgeTest(trainX, trainY) #The ridge regression coefficients under 30 different lambda were obtained for k in range(30): #Traverse all ridge regression coefficients matTestX = np.mat(testX); matTrainX = np.mat(trainX) #Test set meanTrain = np.mean(matTrainX,0) #Test set mean varTrain = np.var(matTrainX,0) #Test set variance matTestX = (matTestX - meanTrain) / varTrain #Test set standardization yEst = matTestX * np.mat(wMat[k,:]).T + np.mean(trainY) #Predict y value according to ws errorMat[i, k] = rssError(yEst.T.A, np.array(testY)) #Statistical error meanErrors = np.mean(errorMat,0) #Calculate the average error of each cross validation minMean = float(min(meanErrors)) #Find the minimum error bestWeights = wMat[np.nonzero(meanErrors == minMean)] #Find the best regression coefficient xMat = np.mat(xArr); yMat = np.mat(yArr).T meanX = np.mean(xMat,0); varX = np.var(xMat,0) unReg = bestWeights / varX #Data is standardized and therefore needs to be restored print('%f%+f*particular year%+f*Number of parts%+f*Is it brand new%+f*original price' % ((-1 * np.sum(np.multiply(meanX,unReg)) + np.mean(yMat)), unReg[0,0], unReg[0,1], unReg[0,2], unReg[0,3])) def ridgeTest(xArr, yArr): """ Function description:Ridge regression test Parameters: xMat - x data set yMat - y data set Returns: wMat - Regression coefficient matrix """ xMat = np.mat(xArr); yMat = np.mat(yArr).T #Data standardization yMean = np.mean(yMat, axis = 0) #Row to row operation, find the mean yMat = yMat - yMean #Data minus mean xMeans = np.mean(xMat, axis = 0) #Row to row operation, find the mean xVar = np.var(xMat, axis = 0) #Row to row operation, variance xMat = (xMat - xMeans) / xVar #Data is normalized by subtracting mean divided by variance numTestPts = 30 #30 different lambda tests wMat = np.zeros((numTestPts, np.shape(xMat)[1])) #Initial regression coefficient matrix for i in range(numTestPts): #Changing lambda to calculate regression coefficient ws = ridgeRegres(xMat, yMat, np.exp(i - 10)) #lambda changes exponentially with e, initially a very small number, wMat[i, :] = ws.T #Calculate regression coefficient matrix return wMat if __name__ == '__main__': lgX = [] lgY = [] setDataCollect(lgX, lgY) crossValidation(lgX, lgY)

The samples are randomly selected here. Because of their randomness, the results of each run may be slightly different. However, as shown in the figure above, it is not very different from the conventional least squares method, that is, ordinary linear regression. We expected to find a more understandable model, which obviously did not achieve the expected effect.

Now let's look at how the regression coefficient changes during the reduction process. The code is as follows:

import numpy as np from bs4 import BeautifulSoup import random def scrapePage(retX, retY, inFile, yr, numPce, origPrc): """ Function description:Read data from the page and generate retX and retY list Parameters: retX - data X retY - data Y inFile - HTML file yr - particular year numPce - Number of LEGO parts origPrc - original price Returns: nothing """ # Open and read HTML file with open(inFile, encoding='utf-8') as f: html = f.read() soup = BeautifulSoup(html) i = 1 # Parse according to HTML page structure currentRow = soup.find_all('table', r = "%d" % i) while(len(currentRow) != 0): currentRow = soup.find_all('table', r = "%d" % i) title = currentRow[0].find_all('a')[1].text lwrTitle = title.lower() # Find out if there are brand new labels if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1): newFlag = 1.0 else: newFlag = 0.0 # To find out whether the logo has been sold, we only collect the data that has been sold soldUnicde = currentRow[0].find_all('td')[3].find_all('span') if len(soldUnicde) == 0: print("commodity #%d not sold "% i) else: # Analyze the page to get the current price soldPrice = currentRow[0].find_all('td')[4] priceStr = soldPrice.text priceStr = priceStr.replace('$','') priceStr = priceStr.replace(',','') if len(soldPrice) > 1: priceStr = priceStr.replace('Free shipping', '') sellingPrice = float(priceStr) # Remove incomplete package price if sellingPrice > origPrc * 0.5: print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice)) retX.append([yr, numPce, newFlag, origPrc]) retY.append(sellingPrice) i += 1 currentRow = soup.find_all('table', r = "%d" % i) def ridgeRegres(xMat, yMat, lam = 0.2): """ Function description:Ridge regression Parameters: xMat - x data set yMat - y data set lam - Reduction factor Returns: ws - regression coefficient """ xTx = xMat.T * xMat denom = xTx + np.eye(np.shape(xMat)[1]) * lam if np.linalg.det(denom) == 0.0: print("The matrix is singular,Cannot transpose") return ws = denom.I * (xMat.T * yMat) return ws def setDataCollect(retX, retY): """ Function description:Read the data of six Lego sets in turn and generate the data matrix Parameters: nothing Returns: nothing """ scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99) #LEGO 8288 in 2006, number of parts 800, original price 49.99 scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99) #LEGO 10030 in 2002, number of parts 3096, original price 269.99 scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99) #LEGO 10179 in 2007, number of parts 5195, original price 499.99 scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99) #LEGO 10181 in 2007, number of parts 3428, original price 199.99 scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99) #LEGO 10189 in 2008, number of parts 5922, original price 299.99 scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99) #LEGO 10196 in 2009, 3263 parts, original price 249.99 def regularize(xMat, yMat): """ Function description:Data standardization Parameters: xMat - x data set yMat - y data set Returns: inxMat - Standardized x data set inyMat - Standardized y data set """ inxMat = xMat.copy() #Data copy inyMat = yMat.copy() yMean = np.mean(yMat, 0) #Row to row operation, find the mean inyMat = yMat - yMean #Data minus mean inMeans = np.mean(inxMat, 0) #Row to row operation, find the mean inVar = np.var(inxMat, 0) #Row to row operation, variance # print(inxMat) print(inMeans) # print(inVar) inxMat = (inxMat - inMeans) / inVar #Data is normalized by subtracting mean divided by variance return inxMat, inyMat def rssError(yArr,yHatArr): """ Function description:Calculate square error Parameters: yArr - Estimate yHatArr - True value Returns: """ return ((yArr-yHatArr)**2).sum() def standRegres(xArr,yArr): """ Function description:Calculate regression coefficient w Parameters: xArr - x data set yArr - y data set Returns: ws - regression coefficient """ xMat = np.mat(xArr); yMat = np.mat(yArr).T xTx = xMat.T * xMat #The regression coefficient is calculated according to the publicity derived in this paper if np.linalg.det(xTx) == 0.0: print("The matrix is singular,Cannot transpose") return ws = xTx.I * (xMat.T*yMat) return ws def ridgeTest(xArr, yArr): """ Function description:Ridge regression test Parameters: xMat - x data set yMat - y data set Returns: wMat - Regression coefficient matrix """ xMat = np.mat(xArr); yMat = np.mat(yArr).T #Data standardization yMean = np.mean(yMat, axis = 0) #Row to row operation, find the mean yMat = yMat - yMean #Data minus mean xMeans = np.mean(xMat, axis = 0) #Row to row operation, find the mean xVar = np.var(xMat, axis = 0) #Row to row operation, variance xMat = (xMat - xMeans) / xVar #Data is normalized by subtracting mean divided by variance numTestPts = 30 #30 different lambda tests wMat = np.zeros((numTestPts, np.shape(xMat)[1])) #Initial regression coefficient matrix for i in range(numTestPts): #Changing lambda to calculate regression coefficient ws = ridgeRegres(xMat, yMat, np.exp(i - 10)) #lambda changes exponentially with e, initially a very small number, wMat[i, :] = ws.T #Calculate regression coefficient matrix return wMat if __name__ == '__main__': lgX = [] lgY = [] setDataCollect(lgX, lgY) print(ridgeTest(lgX, lgY))

Looking at the first line of the run result, you can see that the largest is item 4 and the second largest is item 2.

Therefore, if only one feature is selected for prediction, we should select the fourth feature. If 2 features can be selected, the 4th and 2nd features should be selected.

This analysis method enables us to mine the internal laws of a large number of data. When there are only four features, the effect of this method may not be obvious; But if there are more than 100 features, the method will become very effective: it can point out which features are key and which features are unimportant.

### 6, Summary

- Like classification, regression is also the process of predicting the target value. The difference between regression and classification is that the former predicts continuous variables, while the latter predicts discrete variables.
- Ridge regression is a kind of reduction method, which is equivalent to limiting the size of regression coefficient. Another good reduction method is lasso. Lasso is difficult to solve, but the approximate solution can be obtained by the simple stepwise linear regression method.
- The reduction method can also be regarded as a method to reduce the deviation of a model at the same time.