【翻译自 : How to Use Optimization Algorithms to Manually Fit Regression Models

         【说明:Jason Brownlee PhD大神的文章个人很喜欢,所以闲暇时间里会做一点翻译和学习实践的工作,这里是相应工作的实践记录,希望能帮到有需要的人!】

 

 

        像线性回归和逻辑回归这样的模型都是通过最小二乘优化来训练的,这是找到最小化这些模型误差的系数的最有效方法。但是,可以使用替代的优化算法将回归模型拟合到训练数据集。这对于了解更多有关回归的功能以及应用机器学习中优化的核心性质可能是一个有用的练习。对于不符合最小二乘优化程序要求的数据进行回归,也可能需要它。

        在本教程中,您将发现如何手动优化回归模型的系数。完成本教程后,您将知道:

如何开发推理模型以从头开始进行回归。
如何优化线性回归模型的系数以预测数值。
如何使用随机爬山优化逻辑回归模型的系数。

教程概述

      本教程分为三个部分:他们是:

优化回归模型
优化线性回归模型
优化逻辑回归模型

优化回归模型

        回归模型(如线性回归和逻辑回归)是统计领域公认的算法。两种算法都是线性的,这意味着模型的输出是输入的加权和。线性回归是为需要预测数字的“回归”问题而设计的,逻辑回归是为需要分类标签的“分类”问题而设计的。

这些回归模型涉及使用优化算法来找到模型的每个输入的一组系数,从而将预测误差最小化。由于模型是线性的并且易于理解,因此可以使用有效的优化算法。在线性回归的情况下,可以通过最小二乘法优化来找到系数,这可以使用线性代数来求解。在逻辑回归的情况下,通常使用局部搜索优化算法。可以使用任意优化算法来训练线性和逻辑回归模型。也就是说,我们可以定义一个回归模型,并使用给定的优化算法来找到该模型的一组系数,这些系数会导致最小的预测误差或最大的分类精度。与使用推荐的优化算法相比,使用替代优化算法的平均效率较低。但是,在某些特定情况下,例如输入数据不符合模型的期望(如高斯分布)并且与外部输入不相关时,它可能会更有效。在训练机器学习算法(尤其是回归模型)中展示优化的中心性质也可能是一个有趣的练习。接下来,让我们探讨如何使用随机爬山训练线性回归模型。

优化线性回归模型

          线性回归模型可能是从数据中学习的最简单的预测模型。该模型对每个输入都有一个系数,而预测的输出只是一些输入和系数的权重。

          在本节中,我们将优化线性回归模型的系数。首先,让我们定义一个综合回归问题,我们可以将其用作优化模型的重点。我们可以使用make_regression()函数定义一个具有1,000行和10个输入变量的回归问题。

         下面的示例创建数据集并总结数据的形状。

# define a regression dataset
from sklearn.datasets import make_regression
# define dataset
X, y = make_regression(n_samples=1000, n_features=10, n_informative=2, noise=0.2, random_state=1)
# summarize the shape of the dataset
print(X.shape, y.shape)

       运行示例将打印创建的数据集的形状,从而确认我们的期望。

(1000, 10) (1000,)

        接下来,我们需要定义一个线性回归模型。

       在优化模型系数之前,我们必须开发模型并对其运行方式充满信心。

      首先,我们开发一个函数,该函数针对数据集中给定输入数据行计算模型的激活程度。此函数将获取数据行和模型系数,并通过添加额外的y截距(也称为偏移或偏差)系数来计算输入的加权和。 下面的predict_row()函数实现了这一点。我们使用简单的Python列表和命令式编程风格,而不是NumPy数组或列表压缩,以使代码对于Python初学者更加可读。 随时对其进行优化,并在下面的注释中发布您的代码。

# linear regression
def predict_row(row, coefficients):
	# add the bias, the last coefficient
	result = coefficients[-1]
	# add the weighted input
	for i in range(len(row)):
		result += coefficients[i] * row[i]
	return result

       接下来,我们可以为给定数据集中的每一行调用predict_row()函数。 下面的predict_dataset()函数实现了这一点。

       同样,我们有意使用一种简单的命令式编码方式来提高可读性,而不是使用列表压缩。

# use model coefficients to generate predictions for a dataset of rows
def predict_dataset(X, coefficients):
	yhats = list()
	for row in X:
		# make a prediction
		yhat = predict_row(row, coefficients)
		# store the prediction
		yhats.append(yhat)
	return yhats

       最后,我们可以使用该模型对综合数据集进行预测,以确认其均正常工作。

       我们可以使用rand()函数生成一组随机的模型系数。

       回想一下,我们需要为每个输入(此数据集中的十个输入)增加一个系数,并为y截距系数增加权重。

# define dataset
X, y = make_regression(n_samples=1000, n_features=10, n_informative=2, noise=0.2, random_state=1)
# determine the number of coefficients
n_coeff = X.shape[1] + 1
# generate random coefficients
coefficients = rand(n_coeff)

        然后,我们可以将这些系数与数据集一起使用进行预测。

# generate predictions for dataset
yhat = predict_dataset(X, coefficients)

       我们可以评估这些预测的均方误差。

# calculate model prediction error
score = mean_squared_error(y, yhat)
print('MSE: %f' % score)

      我们可以将所有这些联系在一起,并演示用于回归预测建模的线性回归模型。 下面列出了完整的示例。

# linear regression model
from numpy.random import rand
from sklearn.datasets import make_regression
from sklearn.metrics import mean_squared_error

# linear regression
def predict_row(row, coefficients):
	# add the bias, the last coefficient
	result = coefficients[-1]
	# add the weighted input
	for i in range(len(row)):
		result += coefficients[i] * row[i]
	return result

# use model coefficients to generate predictions for a dataset of rows
def predict_dataset(X, coefficients):
	yhats = list()
	for row in X:
		# make a prediction
		yhat = predict_row(row, coefficients)
		# store the prediction
		yhats.append(yhat)
	return yhats

# define dataset
X, y = make_regression(n_samples=1000, n_features=10, n_informative=2, noise=0.2, random_state=1)
# determine the number of coefficients
n_coeff = X.shape[1] + 1
# generate random coefficients
coefficients = rand(n_coeff)
# generate predictions for dataset
yhat = predict_dataset(X, coefficients)
# calculate model prediction error
score = mean_squared_error(y, yhat)
print('MSE: %f' % score)

        运行示例将为训练数据集中的每个示例生成一个预测,然后为该预测打印均方误差。

       注意:由于算法或评估程序的随机性,或者数值精度的差异,您的结果可能会有所不同。 考虑运行该示例几次并比较平均结果。

       在给定一组随机权重的情况下,我们会期望出现较大的误差,这就是在这种情况下看到的,误差值为大约7,307个单位。

MSE: 7307.756740

       现在,我们可以优化数据集的系数,以在该数据集上实现低误差。

       首先,我们需要将数据集分为训练集和测试集。 重要的是要保留一些在优化模型中不使用的数据,以便在用于对新数据进行预测时,我们可以对模型的性能进行合理的估计。

      我们将67%的数据用于训练,其余33%的数据用作评估模型性能的测试集。

# split into train test sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33)

        接下来,我们可以开发一种随机爬山算法。

        优化算法需要目标函数进行优化。 它必须采用一组系数并返回对应于更好模型的要最小化或最大化的分数。

        在这种情况下,我们将使用给定的系数集评估模型的均方误差,并返回误差分数,必须将其最小化。

        给定数据集和一组系数,下面的Objective()函数将实现此目的,并返回模型的误差。

# objective function
def objective(X, y, coefficients):
	# generate predictions for dataset
	yhat = predict_dataset(X, coefficients)
	# calculate accuracy
	score = mean_squared_error(y, yhat)
	return score

         接下来,我们可以定义随机爬山算法。

         该算法将需要一个初始解(例如随机系数),并将迭代地不断对解做一些小的更改并检查它是否会导致性能更好的模型。 当前解决方案的更改量由step_size超参数控制。 此过程将继续进行固定数量的迭代,也作为超参数提供。

        下面的hillclimbing()函数以数据集,目标函数,初始解和超参数为参数来实现此目的,并返回找到的最佳系数集和估计的性能。

# hill climbing local search algorithm
def hillclimbing(X, y, objective, solution, n_iter, step_size):
	# evaluate the initial point
	solution_eval = objective(X, y, solution)
	# run the hill climb
	for i in range(n_iter):
		# take a step
		candidate = solution + randn(len(solution)) * step_size
		# evaluate candidate point
		candidte_eval = objective(X, y, candidate)
		# check if we should keep the new point
		if candidte_eval <= solution_eval:
			# store the new point
			solution, solution_eval = candidate, candidte_eval
			# report progress
			print('>%d %.5f' % (i, solution_eval))
	return [solution, solution_eval]

        然后我们可以调用此函数,传入一组初始系数作为初始解,并输入训练数据集作为针对模型进行优化的数据集。

# define the total iterations
n_iter = 2000
# define the maximum step size
step_size = 0.15
# determine the number of coefficients
n_coef = X.shape[1] + 1
# define the initial solution
solution = rand(n_coef)
# perform the hill climbing search
coefficients, score = hillclimbing(X_train, y_train, objective, solution, n_iter, step_size)
print('Done!')
print('Coefficients: %s' % coefficients)
print('Train MSE: %f' % (score))

       最后,我们可以评估测试数据集上的最佳模型并报告性能。

# generate predictions for the test dataset
yhat = predict_dataset(X_test, coefficients)
# calculate accuracy
score = mean_squared_error(y_test, yhat)
print('Test MSE: %f' % (score))

       完整实例如下:

# optimize linear regression coefficients for regression dataset
from numpy.random import randn
from numpy.random import rand
from sklearn.datasets import make_regression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

# linear regression
def predict_row(row, coefficients):
	# add the bias, the last coefficient
	result = coefficients[-1]
	# add the weighted input
	for i in range(len(row)):
		result += coefficients[i] * row[i]
	return result

# use model coefficients to generate predictions for a dataset of rows
def predict_dataset(X, coefficients):
	yhats = list()
	for row in X:
		# make a prediction
		yhat = predict_row(row, coefficients)
		# store the prediction
		yhats.append(yhat)
	return yhats

# objective function
def objective(X, y, coefficients):
	# generate predictions for dataset
	yhat = predict_dataset(X, coefficients)
	# calculate accuracy
	score = mean_squared_error(y, yhat)
	return score

# hill climbing local search algorithm
def hillclimbing(X, y, objective, solution, n_iter, step_size):
	# evaluate the initial point
	solution_eval = objective(X, y, solution)
	# run the hill climb
	for i in range(n_iter):
		# take a step
		candidate = solution + randn(len(solution)) * step_size
		# evaluate candidate point
		candidte_eval = objective(X, y, candidate)
		# check if we should keep the new point
		if candidte_eval <= solution_eval:
			# store the new point
			solution, solution_eval = candidate, candidte_eval
			# report progress
			print('>%d %.5f' % (i, solution_eval))
	return [solution, solution_eval]

# define dataset
X, y = make_regression(n_samples=1000, n_features=10, n_informative=2, noise=0.2, random_state=1)
# split into train test sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33)
# define the total iterations
n_iter = 2000
# define the maximum step size
step_size = 0.15
# determine the number of coefficients
n_coef = X.shape[1] + 1
# define the initial solution
solution = rand(n_coef)
# perform the hill climbing search
coefficients, score = hillclimbing(X_train, y_train, objective, solution, n_iter, step_size)
print('Done!')
print('Coefficients: %s' % coefficients)
print('Train MSE: %f' % (score))
# generate predictions for the test dataset
yhat = predict_dataset(X_test, coefficients)
# calculate accuracy
score = mean_squared_error(y_test, yhat)
print('Test MSE: %f' % (score))

       运行该示例将在每次对该模型进行改进时报告迭代次数和均方误差。

       搜索结束时,报告训练数据集上最佳系数集的性能,并计算和报告测试数据集上相同模型的性能。

       注意:由于算法或评估程序的随机性,或者数值精度的差异,您的结果可能会有所不同。 考虑运行该示例几次并比较平均结果。

        在这种情况下,我们可以看到优化算法找到了一组系数,这些系数在训练和测试数据集上均实现了约0.08的误差。

       该算法在训练和测试数据集上发现一个性能非常相似的模型的事实是一个好兆头,表明该模型并未过度拟合(过度优化)训练数据集。 这意味着该模型可以很好地概括新数据。

>1546 0.35426
>1567 0.32863
>1572 0.32322
>1619 0.24890
>1665 0.24800
>1691 0.24162
>1715 0.15893
>1809 0.15337
>1892 0.14656
>1956 0.08042
Done!
Coefficients: [ 1.30559829e-02 -2.58299382e-04  3.33118191e+00  3.20418534e-02
  1.36497902e-01  8.65445367e+01  2.78356715e-02 -8.50901499e-02
  8.90078243e-02  6.15779867e-02 -3.85657793e-02]
Train MSE: 0.080415
Test MSE: 0.080779

         既然我们熟悉了如何手动优化线性回归模型的系数,让我们看一下如何扩展示例以优化用于分类的逻辑回归模型的系数。

优化逻辑回归模型

         Logistic回归模型是线性回归的扩展,用于分类预测模型。

         Logistic回归用于二进制分类任务,这意味着具有两个类标签的数据集为class = 0和class = 1。输出首先涉及计算输入的加权和,然后使该加权和通过逻辑函数(也称为S型函数)。对于属于类别= 1的示例,结果是介于0和1之间的二项式概率。在本节中,我们将基于在上一节中学到的知识来优化用于分类的回归模型的系数。我们将开发模型并使用随机系数对其进行测试,然后使用随机爬坡优化模型系数。首先,让我们定义一个综合二进制分类问题,我们可以将其用作优化模型的重点。我们可以使用make_classification()函数来定义一个包含1,000行和五个输入变量的二进制分类问题。

       下面的示例创建数据集并总结数据的形状。

# define a binary classification dataset
from sklearn.datasets import make_classification
# define dataset
X, y = make_classification(n_samples=1000, n_features=5, n_informative=2, n_redundant=1, random_state=1)
# summarize the shape of the dataset
print(X.shape, y.shape)

         运行示例将打印创建的数据集的形状,从而确认我们的期望。

(1000, 5) (1000,)

        接下来,我们需要定义一个逻辑回归模型。

        让我们首先更新predict_row()函数,以通过逻辑函数传递输入和系数的加权和。

函数功能定义为:

                              logistic = 1.0 /(1.0 + exp(-result))
         其中,结果是输入的加权总和,系数和exp()是e(欧拉数),提高到通过exp()函数实现的提供值的幂。

        下面列出了更新的predict_row()函数。

# logistic regression
def predict_row(row, coefficients):
	# add the bias, the last coefficient
	result = coefficients[-1]
	# add the weighted input
	for i in range(len(row)):
		result += coefficients[i] * row[i]
	# logistic function
	logistic = 1.0 / (1.0 + exp(-result))
	return logistic

        从线性回归到逻辑回归的变化就可以了。

        与线性回归一样,我们可以使用一组随机模型系数来测试模型。

# determine the number of coefficients
n_coeff = X.shape[1] + 1
# generate random coefficients
coefficients = rand(n_coeff)
# generate predictions for dataset
yhat = predict_dataset(X, coefficients)

       该模型做出的预测是属于类别= 1的示例的概率。

      对于预期的类别标签,我们可以将预测取整为整数值0和1。

# round predictions to labels
yhat = [round(y) for y in yhat]

       我们可以评估这些预测的分类准确性。

# calculate accuracy
score = accuracy_score(y, yhat)
print('Accuracy: %f' % score)

        我们可以将所有这些结合在一起,并演示用于二进制分类的简单逻辑回归模型。 下面列出了完整的示例。

# logistic regression function for binary classification
from math import exp
from numpy.random import rand
from sklearn.datasets import make_classification
from sklearn.metrics import accuracy_score

# logistic regression
def predict_row(row, coefficients):
	# add the bias, the last coefficient
	result = coefficients[-1]
	# add the weighted input
	for i in range(len(row)):
		result += coefficients[i] * row[i]
	# logistic function
	logistic = 1.0 / (1.0 + exp(-result))
	return logistic

# use model coefficients to generate predictions for a dataset of rows
def predict_dataset(X, coefficients):
	yhats = list()
	for row in X:
		# make a prediction
		yhat = predict_row(row, coefficients)
		# store the prediction
		yhats.append(yhat)
	return yhats

# define dataset
X, y = make_classification(n_samples=1000, n_features=5, n_informative=2, n_redundant=1, random_state=1)
# determine the number of coefficients
n_coeff = X.shape[1] + 1
# generate random coefficients
coefficients = rand(n_coeff)
# generate predictions for dataset
yhat = predict_dataset(X, coefficients)
# round predictions to labels
yhat = [round(y) for y in yhat]
# calculate accuracy
score = accuracy_score(y, yhat)
print('Accuracy: %f' % score)

       运行示例将为训练数据集中的每个示例生成一个预测,然后打印该预测的分类准确性。

       注意:由于算法或评估程序的随机性,或者数值精度的差异,您的结果可能会有所不同。 考虑运行该示例几次并比较平均结果。

        在给定一组随机权重和每个类中具有相同数量示例的数据集的情况下,我们可以期望达到约50%的准确性,在这种情况下,这大约就是我们看到的结果。

Accuracy: 0.540000

        现在,我们可以优化数据集的权重,以在该数据集上实现良好的准确性。

        用于线性回归的随机爬山算法可以再次用于逻辑回归。

       重要的区别是对objective()函数的更新,以使预测变圆并使用分类精度而不是均方误差来评估模型。

# objective function
def objective(X, y, coefficients):
	# generate predictions for dataset
	yhat = predict_dataset(X, coefficients)
	# round predictions to labels
	yhat = [round(y) for y in yhat]
	# calculate accuracy
	score = accuracy_score(y, yhat)
	return score

         还必须更新hillclimbing()函数以最大化解决方案的得分,而不是在线性回归的情况下最小化。

# hill climbing local search algorithm
def hillclimbing(X, y, objective, solution, n_iter, step_size):
	# evaluate the initial point
	solution_eval = objective(X, y, solution)
	# run the hill climb
	for i in range(n_iter):
		# take a step
		candidate = solution + randn(len(solution)) * step_size
		# evaluate candidate point
		candidte_eval = objective(X, y, candidate)
		# check if we should keep the new point
		if candidte_eval >= solution_eval:
			# store the new point
			solution, solution_eval = candidate, candidte_eval
			# report progress
			print('>%d %.5f' % (i, solution_eval))
	return [solution, solution_eval]

        最后,可以在运行结束时使用分类准确性评估通过搜索找到的系数。

# generate predictions for the test dataset
yhat = predict_dataset(X_test, coefficients)
# round predictions to labels
yhat = [round(y) for y in yhat]
# calculate accuracy
score = accuracy_score(y_test, yhat)
print('Test Accuracy: %f' % (score))

       综上所述,下面列出了使用随机爬坡最大化logistic回归模型的分类准确性的完整示例。

# optimize logistic regression model with a stochastic hill climber
from math import exp
from numpy.random import randn
from numpy.random import rand
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

# logistic regression
def predict_row(row, coefficients):
	# add the bias, the last coefficient
	result = coefficients[-1]
	# add the weighted input
	for i in range(len(row)):
		result += coefficients[i] * row[i]
	# logistic function
	logistic = 1.0 / (1.0 + exp(-result))
	return logistic

# use model coefficients to generate predictions for a dataset of rows
def predict_dataset(X, coefficients):
	yhats = list()
	for row in X:
		# make a prediction
		yhat = predict_row(row, coefficients)
		# store the prediction
		yhats.append(yhat)
	return yhats

# objective function
def objective(X, y, coefficients):
	# generate predictions for dataset
	yhat = predict_dataset(X, coefficients)
	# round predictions to labels
	yhat = [round(y) for y in yhat]
	# calculate accuracy
	score = accuracy_score(y, yhat)
	return score

# hill climbing local search algorithm
def hillclimbing(X, y, objective, solution, n_iter, step_size):
	# evaluate the initial point
	solution_eval = objective(X, y, solution)
	# run the hill climb
	for i in range(n_iter):
		# take a step
		candidate = solution + randn(len(solution)) * step_size
		# evaluate candidate point
		candidte_eval = objective(X, y, candidate)
		# check if we should keep the new point
		if candidte_eval >= solution_eval:
			# store the new point
			solution, solution_eval = candidate, candidte_eval
			# report progress
			print('>%d %.5f' % (i, solution_eval))
	return [solution, solution_eval]

# define dataset
X, y = make_classification(n_samples=1000, n_features=5, n_informative=2, n_redundant=1, random_state=1)
# split into train test sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33)
# define the total iterations
n_iter = 2000
# define the maximum step size
step_size = 0.1
# determine the number of coefficients
n_coef = X.shape[1] + 1
# define the initial solution
solution = rand(n_coef)
# perform the hill climbing search
coefficients, score = hillclimbing(X_train, y_train, objective, solution, n_iter, step_size)
print('Done!')
print('Coefficients: %s' % coefficients)
print('Train Accuracy: %f' % (score))
# generate predictions for the test dataset
yhat = predict_dataset(X_test, coefficients)
# round predictions to labels
yhat = [round(y) for y in yhat]
# calculate accuracy
score = accuracy_score(y_test, yhat)
print('Test Accuracy: %f' % (score))

         每次对模型进行改进时,运行示例将报告迭代次数和分类准确性。

         搜索结束时,报告训练数据集上最佳系数集的性能,并计算和报告测试数据集上相同模型的性能。

        注意:由于算法或评估程序的随机性,或者数值精度的差异,您的结果可能会有所不同。 考虑运行该示例几次并比较平均结果。

        在这种情况下,我们可以看到优化算法找到了一组权重,这些权重在训练数据集上达到约87.3%的准确度,在测试数据集上达到约83.9%的准确度。

>200 0.85672
>225 0.85672
>230 0.85672
>245 0.86418
>281 0.86418
>285 0.86716
>294 0.86716
>306 0.86716
>316 0.86716
>317 0.86716
>320 0.86866
>348 0.86866
>362 0.87313
>784 0.87313
>1649 0.87313
Done!
Coefficients: [-0.04652756  0.23243427  2.58587637 -0.45528253 -0.4954355  -0.42658053]
Train Accuracy: 0.873134
Test Accuracy: 0.839394