概率统计与随机过程--作业5

一、推导题

二、计算题

1、某单位为了研究太阳镜销售和广告费用之间的关系，搜集了以下数据，使用回归分析方法得到线性回归模型：

解：

（1）绘制的散点图和回归线如下图所示：

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']  # 正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

#数据
data=np.array([[2,5,6,7,22,25,28,30,22,18],[75, 90,148,183,242,263,278,318,256,200]])

x=data[0]
y=data[1]
plt.scatter(x, y,  c='r',marker='o',label='销售量') #散点图
linreg = LinearRegression()#线性回归
linreg.fit(x.reshape(-1,1),y) #拟合,x要转换为列向量
y_pre=linreg.predict(x.reshape(-1,1))
plt.plot(x,y_pre,c='b') #回归线
s='y='
for i in range(len(linreg.coef_)):
    if(linreg.coef_[i]>=0 and i>0):
        s=s+'+'+str(round(linreg.coef_[i],3))+'x'+str(i)
    else:
        s=s+str(round(linreg.coef_[i],3))+'x'+str(i)
if(linreg.intercept_>=0):
    s=s+'+'+str(round(linreg.intercept_,3))
else:
        s=s+str(round(linreg.intercept_,3))
plt.title('太阳镜销售和广告费用之间的关系')
plt.xlabel('x-广告费用（万元）')
plt.ylabel('y-销售量（个）')
plt.legend()
plt.show()
a_i=linreg.intercept_  # a 的估计值
b_i=linreg.coef_[0]    # a 的估计值
print("线性回归方程为：",s)

# 计算统计量
Sxx=0
Syy=0
Sxy=0
SSe=Qe=0

n=data.shape[1]
x_=x.mean()
y_=y.mean()
for i in range(n):
    t=(x[i]-x_)**2
    Sxx=Sxx+t
    t=(y[i]-y_)**2
    Syy=Syy+t
    t=(y[i]-y_)*(x[i]-x_)
    Sxy=Sxy+t
    t=(y[i]-y_pre[i])**2
    SSe=SSe+t
# b的估计值 b_i=Sxy/Sxx
Qe=Syy-b_i*Sxy # Qe==SSe
sigma_i= np.sqrt(Qe/(n-2))  #sigma 的估计值
print("主要统计参数：Sxx={:.3f},Syy={:.3f},Sxy={:.3f},SSe=Qe={:.3f},Sigma={:.3f}".format(Sxx,Syy,Sxy,SSe,sigma_i))
sigma_i= np.sqrt(Qe/(n-2))  #sigma 的估计值
x_i=35 #输入的x值
y_i=b_i*x_i+a_i #相应的预测值
t_c=2.306 # t_a/2的临界值，a=0.05
interval=np.sqrt((1+1/n+(x_i-x_)**2/Sxx)*sigma**2)*t_c
print("对应x={:.3f}的Y的预测值为{:.3f}，置信度为95%的预测区间为:({:.3f}，{:.3f})".format(x_i,y_i,y_i-interval,y_i+interval))

 2. 对鲍鱼数据集（abalone.txt）进行向前逐步回归，将“Length”列值全设置为1，给出优化后属性列表（参加ppt中【例5-9】，【例5-10】及相关代码）。

答案： final formula is Age~Rings+Viscera+Height+Shucked+Shell+Whole

import numpy as np
import pandas as pd
import statsmodels.api as sm #最小二乘
from statsmodels.formula.api import ols #加载ols模型
# 数据准备
#读取鲍鱼数据集
aba = pd.read_table('abalone.txt',sep=',', names=['Length',	'Diam',	'Height',	'Whole'	,'Shucked',	'Viscera',	'Shell',	'Rings','Age'
],header = None)#该数据集源于UCI，记录了鲍⻥的⽣物属性，⽬标字段是该⽣物的年龄
print(aba.shape)
aba.iloc[:, 0] = 1 # 把类型列置1
print(aba.head())

print(aba.shape) #查看数据集大小
print(aba.head(5)) #查看前10行数据
print(aba.columns)

#定义向前逐步回归函数
def forward_select(data,target):
    variate=set(data.columns)  #将字段名转换成字典类型
    variate.remove(target)  #去掉因变量的字段名
    selected=[]
    current_score,best_new_score=float('inf'),float('inf')  #目前的分数和最好分数初始值都为无穷大（因为AIC越小越好）
    #循环筛选变量
    while variate:
        aic_with_variate=[]
        for candidate in variate:  #逐个遍历自变量
            formula="{}~{}".format(target,"+".join(selected+[candidate]))  #将自变量名连接起来
            aic=ols(formula=formula,data=data).fit().aic  #利用ols训练模型得出aic值
            aic_with_variate.append((aic,candidate))  #将第每一次的aic值放进空列表
        aic_with_variate.sort(reverse=True)  #降序排序aic值
        best_new_score,best_candidate=aic_with_variate.pop()  #最好的aic值等于删除列表的最后一个值，以及最好的自变量等于列表最后一个自变量
        if current_score>best_new_score:  #如果目前的aic值大于最好的aic值
            variate.remove(best_candidate)  #移除加进来的变量名，即第二次循环时，不考虑此自变量了
            selected.append(best_candidate)  #将此自变量作为加进模型中的自变量
            current_score=best_new_score  #最新的分数等于最好的分数
            print("aic is {},continuing!".format(current_score))  #输出最小的aic值
        else:
            print("for selection over!")
            break
    formula="{}~{}".format(target,"+".join(selected))  #最终的模型式子
    print("final formula is {}".format(formula))
    model=ols(formula=formula,data=data).fit()
    return(model)
# 对数据进行前向逐步回归
forward_select(data=aba,target="Age")