最小二乘回归stata命令

房屋面积和售价

线性回归：首先一个例子是房屋面积和价格，这个就构成了一个简单的线性回归

import pandas as pd
import numpy as np
from io import StringIO
from sklearn import linear_model
import matplotlib.pyplot as plt
csv_data = 'square_feet,price\n150,6450\n200,7450\n250,8450\n300,9450\n350,11450\n400,15450\n600,18450\n'
#<class 'pandas.core.frame.DataFrame'>这个是pandas里面的共用的数据类型
df = pd.read_csv(StringIO(csv_data))
print(type(df))
#创建实例
regr = linear_model.LinearRegression()
#x输入面积，y是价格
#进行拟合
#df['square_feet'].values.reshape(-1,1)这个是输入参数98
regr.fit(df['square_feet'].values.reshape(-1,1), df['price'])
a, b = regr.coef_, regr.intercept_
print(a)
#截距表示什么呢？就是b啊
print(b)
area = 238.5
print(a*area+b)
#输入一个值，拿来预测这个结果
regr.predict(area)
#这个plt的scatter函数来绘图,这个是绘制散列点
plt.scatter(df['square_feet'], df['price'])
#这个是画线
plt.plot(df['square_feet'], regr.predict(df['square_feet'].values.reshape(-1,1)))
#plt.show()
#预测的y的数组
regr.predict(df['square_feet'].values.reshape(-1,1))

结果：

如何拟合计算出m和b

线性回归的标准是所有实际的点到这个直线的距离最短：

将绝对值计算成平方，然后求和如下：

这个就是常说的损失函数，Xi和Yi还有N都是已知数，那么未知就是m和b，实际转换一下就是求一个二元二次方程的最小值。

补充：

Z = X*X+Y*Y+2*Y

绘制图像如下：

from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from matplotlib import pyplot as plt
 
fig = plt.figure()
ax = Axes3D(fig)
x=np.arange(-5,5,0.1)
y=np.arange(-5,5,0.1)
X, Y = np.meshgrid(x, y)#网格的创建，这个是关键
Z=X*X+Y*Y+2*Y
plt.xlabel('x')
plt.ylabel('y')
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow')
plt.show()

单从一个面上来看是不是平的点就是最小点呢！(这个主要是为下面打基础)

最小二乘法

一般最小化误差的计算有两种方法：最小二乘法+梯度下降法

对m和b求偏导，然后令他等于0

转换成二元一次方程求解，这样就可以直接计算出m和b，就得到那个直线了。

梯度下降法：

微分：微分有两种说法，某点的切线的斜率/函数的变化率

梯度：梯度实际上就是多变量微分的一般化，梯度就是分别对每个变量进行微分，然后用逗号分割开，梯度是用<>包括起来，说明梯度其实一个向量。

• 在单变量的函数中，梯度其实就是函数的微分，代表着函数在某个给定点的切线的斜率
• 在多变量函数中，梯度是一个向量，向量有方向，梯度的方向就指出了函数在给定点的上升最快的方向

梯度下降的数学公式：

J是关于Θ的一个函数，我们当前所处的位置为Θ0点，要从这个点走到J的最小值点，也就是山底。首先我们先确定前进的方向，也就是梯度的反向，然后走一段距离的步长，也就是α，走完这个段步长，就到达了Θ1这个点

α：学习率或者叫做步长

import numpy as np
import pylab

def compute_error(b,m,data):

    totalError = 0
    #Two ways to implement this
    #first way
    # for i in range(0,len(data)):
    #     x = data[i,0]
    #     y = data[i,1]
    #
    #     totalError += (y-(m*x+b))**2

    #second way
    x = data[:,0]
    y = data[:,1]
    totalError = (y-m*x-b)**2
    totalError = np.sum(totalError,axis=0)

    return totalError/float(len(data))

def optimizer(data,starting_b,starting_m,learning_rate,num_iter):
    b = starting_b
    m = starting_m

    #gradient descent
    for i in range(num_iter):
        #update b and m with the new more accurate b and m by performing
        # thie gradient step
        b,m =compute_gradient(b,m,data,learning_rate)
        if i%100==0:
            print('iter {0}:error={1}'.format(i,compute_error(b,m,data)))
    return [b,m]

def compute_gradient(b_current,m_current,data ,learning_rate):

    b_gradient = 0
    m_gradient = 0

    N = float(len(data))
    #Two ways to implement this
    #first way
    #print(N)
    for i in range(0,len(data)):
        x = data[i,0]
        y = data[i,1]
    
        #computing partial derivations of our error function
        #b_gradient = -(2/N)*sum((y-(m*x+b))^2)
        #m_gradient = -(2/N)*sum(x*(y-(m*x+b))^2)
        b_gradient += -(2/N)*(y-((m_current*x)+b_current))
        m_gradient += -(2/N) * x * (y-((m_current*x)+b_current))
    
    """
    #Vectorization implementation
    x = data[:,0]
    y = data[:,1]
    b_gradient = -(2/N)*(y-m_current*x-b_current)
    b_gradient = np.sum(b_gradient,axis=0)
    m_gradient = -(2/N)*x*(y-m_current*x-b_current)
    m_gradient = np.sum(m_gradient,axis=0)
        #update our b and m values using out partial derivations
    """
    new_b = b_current - (learning_rate * b_gradient)
    new_m = m_current - (learning_rate * m_gradient)
    return [new_b,new_m]


def plot_data(data,b,m):

    #plottting
    x = data[:,0]
    y = data[:,1]
    y_predict = m*x+b
    pylab.plot(x,y,'o')
    pylab.plot(x,y_predict,'k-')
    pylab.show()


def Linear_regression():
    # get train data
    data =np.loadtxt('data.csv',delimiter=',')

    #define hyperparamters
    #learning_rate is used for update gradient
    #defint the number that will iteration
    # define  y =mx+b
    learning_rate = 0.001
    initial_b =0.0
    initial_m = 0.0
    num_iter = 1000

    #train model
    #print b m error
    print('initial variables:\n initial_b = {0}\n intial_m = {1}\n error of begin = {2} \n'\
    .format(initial_b,initial_m,compute_error(initial_b,initial_m,data)))

    #optimizing b and m
    [b ,m] = optimizer(data,initial_b,initial_m,learning_rate,num_iter)

    #print final b m error
    print('final formula parmaters:\n b = {1}\n m={2}\n error of end = {3} \n'.format(num_iter,b,m,compute_error(b,m,data)))

    #plot result
    plot_data(data,b,m)

if __name__ =='__main__':
    Linear_regression()

运行结果

Error实际上是计算的损失函数的结果