概率统计模型

更新于 分类于 阅读次数： 本文字数： 3.2k 阅读时长 ≈ 5 分钟
用最小二乘法等进行拟合

第七章-概率统计模型

习题1

\(\mu\)的一个置信水平为\(1-\alpha\)的置信区间为

显著水平为\(\alpha =0.1\),\(\alpha /2=0.05\),\(n=5\)，服从正态分布，则其由给出的数据算得均值和方差，最后计算出\(\mu\)的置信水平为0.90的区间。

代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
d = np.array([1050,1100,1120,1250,1280])
n = len(d)
# 计算均值
xb = d.mean()
# 计算标准差
s = d.std(ddof=1)
# 计算样本均值的标准差
sm = sem(d)
a = 0.10
ta = t.ppf(1-a/2,n-1)
L = [xb-sm*ta,xb+sm*ta]
print(L)

结果

1
L=[1064.89955998421, 1255.10044001579]

习题2-1

构建样本点

储存在points矩阵中。

1
2
3
4
5
6
7
8
9
10
11
x=np.linspace(0,10,N)
a = 1.2
b = 3.5
y0 = a * x + b
y_noise = [np.random.normal(0, 0.1) + y for y in y0]
points = []
i = 0
for x0 in x :
    points.append([x0,np.random.normal(0, 0.1) + a * x0 + b])
    i = i + 1
points = np.array(points)

leastsq拟合

先尝试用sklearn的leastsq拟合（此方法用的是最小二乘）得到结果作为验证。

1
2
3
4
5
p0 = np.random.randn(2)
# t=[a,b]
fun = lambda t, x : t[0] * x + t[1]
err = lambda t, x, z: fun(t,x) - z
p2 = leastsq(err, p0, args=(x,y_noise))[0]

结果：

1
2
3
leastsq拟合: [1.20330845 3.48951605]
a = 1.20330845
b = 3.48951605

梯度下降法

接下去分别用设置不同的损失函数实现最小二乘法、最小绝对偏离估计、最小最大偏离估计用于线性回归的分析。

定义损失函数

最小二乘法

对应的损失函数为：

\[ f=\sum (ax_j+b_j-y_j)^2 \]

1
2
3
4
5
6
7
8
9
10
def computer_cost(w, b, points):
    total_cost = 0
    M = len(points)
    # 逐点计算平方损失误差，然后求平均数
    for i in range(M):
        x = points[i][0]
        y = points[i][1]
        total_cost += (y - w * x - b) ** 2
    # 取平均
    return total_cost / M

最小绝对偏离估计

损失函数：

\[ f=\sum|sx_j+b_j-y_j| \]

1
2
3
4
5
6
7
8
9
10
def computer_cost_B(w, b, points):
    total_cost = 0
    M = len(points)
    # 逐点计算平方损失误差，然后求平均数
    for i in range(M):
        x = points[i][0]
        y = points[i][1]
        total_cost += abs(y - w * x - b)
    # 取平均
    return total_cost / M

最小最大偏离估计

损失函数：

\[ f=max|sx_j+b_j-y_j| \]

1
2
3
4
5
6
7
8
9
10
11
12
13
def computer_cost_C(w, b, points):
    total_cost = 0
    M = len(points)
    # 逐点计算平方损失误差，然后求平均数
    maxCost = 0
    for i in range(M):
        x = points[i][0]
        y = points[i][1]
        if abs(y - w * x - b) > maxCost:
            maxCost = abs(y - w * x - b)
        else:
            continue
    return maxCost

定义迭代参数

1
2
3
4
5
6
7
8
# 步长
alpha = 0.0001
# 初始w
initial_w = 1
# 初始b
initial_b = 1
# 迭代次数
num_iter = 1000

由于迭代下降最终的a和b结果的迭代结果在不同的损失函数下都是相同的，所以当cost小于等于1的时候进行截断，观察不同方法的拟合效果：

可见当前迭代未结束，相同cost情况下最小最大偏离估计的拟合速度更快：

1
2
3
4
5
6
7
8
9
10
11
12
13
leastsq拟合: [1.20330845 3.48951605]
最小二乘法:
a is : 1.4119904872895737
b is : 2.1132899832906866
cost is : 0.49050095850618497
最小绝对偏离估计:
a is : 1.3806101062114216
b is : 2.321991385922966
cost is : 0.49535272270004677
最小最大偏离估计:
a is : 1.2607599489963919
b is : 3.119078452972342
cost is : 0.4972349489227139

结果

1
2
3
4
5
6
7
8
9
10
11
12
13
leastsq拟合: [1.20330845 3.48951605]
最小二乘法:
a is : 1.207278913671576
b is : 3.4747646079211276
cost is : 0.008561611307963823
最小绝对偏离估计:
a is : 1.207278913671576
b is : 3.4747646079211276
cost is : 0.07193785933645268
最小最大偏离估计:
a is : 1.207278913671576
b is : 3.4747646079211276
cost is : 0.25574581839753696

可见最小二乘法拟合效果最好。

最终拟合结果，三种方法近乎相同：

习题2-2

随机将百分之十的数据y加减20后。

在cost为50的时候进行截断，可以发现第二种方法的拟合速度较快。

最终的拟合结果：

1
2
3
4
5
6
7
8
# 步长
alpha = 0.001
# 初始w
initial_w = 1
# 初始b
initial_b = 3
# 迭代次数
num_iter = 1000
1
2
3
4
5
6
7
8
9
10
11
12
13
leastsq拟合: [1.35707397 3.52990959]
最小二乘法:
a is : 1.3728617616388785
b is : 3.424909795319972
cost is : 39.02051950964002
最小绝对偏离估计:
a is : 1.3728617616388785
b is : 3.424909795319972
cost is : 2.647783487811292
最小最大偏离估计:
a is : 1.3728617616388785
b is : 3.424909795319972
cost is : 21.373929585192865

可见最小绝对偏离估计的cost最小，拟合速度最快，拟合效果也最好。但是另外两个函数对应的拟合损失较大。

改进方法

最小二乘法的方差能达到很小，但其估计量的方差不一定最小，于是是否可以找到一个有偏估计，使得虽然有微小的偏差，但其精度却能够大大高于无偏的估计量。其实就是岭回归，即将最小二乘法加一个范数，并乘以一定的权重，但能比最小二乘法稳定很多。

1
2
3
4
5
6
7
8
9
dataMat = np.array(points)
X = dataMat[:,0:1]   # 变量x
y = dataMat[:,1]   #变量y
model = Ridge(alpha=0.2)
model = RidgeCV(alphas=[0.1, 1.0, 500])  # 通过RidgeCV可以设置多个参数值，算法使用交叉验证获取最佳参数值
model.fit(X, y)   # 线性回归建模
print('系数矩阵:\n',model.coef_)
print('线性回归模型:\n',model)

其结果如下图所示，可见可以达到很好的拟合效果：