4.1 线性模型的参数估计

4.1线性模型的参数估计案例引入工程师要研究尖卵形弹头的曲径比与钻地弹侵彻深度之间的关系。考虑曲径比这个因子的4个水平：2.5、3、3.5、4，每个水平重复5次侵彻试验得到数据如下：回顾：在第二章中，我们如何分析上述数据？例：弹头曲径比对侵彻深度的影响曲径比观测值总和平均值123452.55.755.425.305.395.7027.565.5123.05.655.935.905.796.1029.375.8743.56.006.516.106.376.2931.276.2544.07.257.007.156.857.1035.357.070

案例引入侵彻深度与曲径比之间可能存在一次或二次多项式的关系，如果能够得到这种潜在规律，就可以实现预测与优化！案例引入某型预制破片弹主要用于毁伤敌方有生力量和轻型装甲目标。为分析炸高与落角对该型预制破片弹的威力的影响，选定8枚该型弹开展地面静态试验。将破片弹固定在设定的高度上以设定的弹着角引爆，测量其对地面的毁伤面积，数据如下：回顾：第三章中，我们如何分析上述数据？例：预制破片弹毁伤面积处理记号试验方案两次试验结果之和122m45°2m90°4m45°4m90°

案例引入固定效应模型主要用于分析定性因子对连续响应的影响;对于定量因子和连续响应,固定效应无法有效解决预测、优化、试验数据丢失等问题!回归模型：定量因子与连续响应之间的函数依赖关系，可用于实现预测、优化等目的。钻地弹侵彻深度数据的再分析预制破片弹毁伤数据的再分析总结提升：线性模型参数估计知识目标：能说清线性回归模型“参数线性”的核心本质；能讲透最小二乘法原理，掌握其矩阵求解步骤。能力目标：面对一元/二元连续因子数据，能独立构建回归模型；能手动或用Python完成参数估计求解；能解读参数估计结果，并初步分析模型适配性。学习目标钻地弹侵彻深度数据的再分析预制破片弹毁伤数据的再分析总结提升：线性模型参数估计一.侵彻深度数据的再分析首先需要“猜测”一个合适的模型。一种方法是依靠物理直觉构建曲径比与侵彻深度之间的函数关系。另一种方法是从散点图中获取灵感。（一）模型构建一.侵彻深度数据的再分析

（一）模型构建一.侵彻深度数据的再分析

（二）最小二乘估计探究一.侵彻深度数据的再分析

（二）最小二乘估计探究

如何求解？一.侵彻深度数据的再分析已知

（二）最小二乘估计探究一.侵彻深度数据的再分析importnumpy

asnp

importmatplotlib.pyplot

asplt

#1.输入原始数据

x_levels=[2.5,3.0,3.5,4.0]

y_data=[

[5.75,5.42,5.30,5.39,5.70],#x=2.5

[5.65,5.93,5.90,5.79,6.10],#x=3.0

[6.00,6.51,6.10,6.37,6.29],#x=3.5

[7.25,7.00,7.15,6.85,7.10]#x=4.0

]

x_list=[]

y_list=[]

repeats=5

fori,xinenumerate(x_levels):

forjinrange(repeats):

x_list.append(x)

y_list.append(y_data[i][j])

x=np.array(x_list)#(20,)

y=np.array(y_list)#(20,)#2.构造设计矩阵X:[1,x,x^2]

X=np.column_stack([np.ones(len(x)),x,x**2])

beta_hat=np.linalg.inv(X.T@X)@X.T@y

beta_0,beta_1,beta_2=beta_hat

#3.计算RSS和误差方差估计

y_fitted=X@beta_hat

residuals=y-y_fitted

RSS=np.sum(residuals**2)

n=len(y)

p=X.shape[1]

df_error=n-p

MSE=RSS/df_error

#输出结果（可选打印）

print(f"回归系数估计：β₀={beta_0:.4f},β₁={beta_1:.4f},β₂={beta_2:.4f}")

print(f"残差平方和(RSS)={RSS:.6f}")

print(f"误差方差估计σ²_hat={MSE:.6f}")

（二）最小二乘估计探究一.侵彻深度数据的再分析

你能预测吗？预测准确吗？哪些因素会影响预测准确性？一.侵彻深度数据的再分析模型预测效果取决于两个方面：认知不确定性：模型本身是否合适，即是否有模型误差的问题；随机不确定性：每次试验都带随机性，需要量化这种随机性！

（三）模型评估与改进一.侵彻深度数据的再分析（三）模型评估与改进一.侵彻深度数据的再分析

（三）模型评估与改进一.侵彻深度数据的再分析importnumpy

asnp

importmatplotlib.pyplot

asplt

#1.输入原始数据

x_levels=[2.5,3.0,3.5,4.0]

y_data=[

[5.75,5.42,5.30,5.39,5.70],#x=2.5

[5.65,5.93,5.90,5.79,6.10],#x=3.0

[6.00,6.51,6.10,6.37,6.29],#x=3.5

[7.25,7.00,7.15,6.85,7.10]#x=4.0

]

x_list=[]

y_list=[]

repeats=5

fori,xinenumerate(x_levels):

forjinrange(repeats):

x_list.append(x)

y_list.append(y_data[i][j])

x=np.array(x_list)#(20,)

y=np.array(y_list)#(20,)#2.求回归系数的最小二乘估计X=np.column_stack([np.ones(len(x)),x])

beta_hat=np.linalg.inv(X.T@X)@X.T@y

如何修改代码？X=np.column_stack([np.ones(len(x)),x,

x**2])

（三）模型评估与改进一.侵彻深度数据的再分析（三）模型评估与改进一.侵彻深度数据的再分析

试验次数模型参数个数试验次数减去线性回归模型中参数的个数为残差平方和的自由度！（四）误差方差的估计一.侵彻深度数据的再分析importnumpy

asnp

importmatplotlib.pyplot

asplt

#1.输入原始数据

x_levels=[2.5,3.0,3.5,4.0]

y_data=[

[5.75,5.42,5.30,5.39,5.70],#x=2.5

[5.65,5.93,5.90,5.79,6.10],#x=3.0

[6.00,6.51,6.10,6.37,6.29],#x=3.5

[7.25,7.00,7.15,6.85,7.10]#x=4.0

]

x_list=[]

y_list=[]

repeats=5

fori,xinenumerate(x_levels):

forjinrange(repeats):

x_list.append(x)

y_list.append(y_data[i][j])

x=np.array(x_list)#(20,)

y=np.array(y_list)#(20,)

#2.求回归系数最小二乘估计

X=np.column_stack([np.ones(len(x)),x**2])

beta_hat=np.linalg.inv(X.T@X)@X.T@y（四）误差方差的估计一.侵彻深度数据的再分析

（四）误差方差的估计钻地弹侵彻深度数据的再分析预制破片弹毁伤数据的再分析总结提升：线性模型参数估计案例引入某型预制破片弹主要用于毁伤敌方有生力量和轻型装甲目标。为分析炸高与落角对该型预制破片弹的威力的影响，选定8枚该型弹开展地面静态试验。将破片弹固定在设定的高度上以设定的弹着角引爆，测量其对地面的毁伤面积，数据如下：回顾：第三章中，我们如何分析上述数据？例：预制破片弹毁伤面积试验方案两次试验结果之和122m45°2m90°4m45°4m90°

二.预制破片弹毁伤数据的再分析试验方案两次试验结果之和122m45°2m90°4m45°4m90°炸高和落角都是定量因子，采用回归模型对数据进行建模。（一）模型构建二.预制破片弹毁伤数据的再分析试验方案编码变量两次试验结果之和122m45°2m90°4m45°4m90°

编码变换有何好处？构建一个什么样的回归模型？（一）模型构建二.预制破片弹毁伤数据的再分析

（一）模型构建

试验方案编码变量两次试验结果之和122m45°2m90°4m45°4m90°二.预制破片弹毁伤数据的再分析

（二）模型参数估计

二.预制破片弹毁伤数据的再分析

你能解释这些模型参数的含义吗？

（二）模型参数估计

二.预制破片弹毁伤数据的再分析

例：预制破片弹毁伤面积处理记号试验方案两次试验结果之和122m45°2m90°4m45°4m90°二.预制破片弹毁伤数据的再分析问题一：若第7号试验数据丢失，还能拟合回归模型吗？importnumpyasnp

data=np.array([

[-1,-1,1202.44],

[-1,-1,1509.92],

[-1,1,3171.77],

[-1,1,3486.43],

[1,-1,2742.72],

[1,-1,2423.20],

[1,1,2501.86],

[1,1,3155.33]])

x1_all=data[:,0]

x2_all=data[:,1]

y_all=data[:,2]

X=np.column_stack([np.ones(len(y_all),x1_all,x2_all,

x1_all*x2_all])

#最小二乘估计:β=(X^TX)^{-1}X^Ty

beta_hat=np.linalg.inv(X.T@X)@X.T@y_all

二.预制破片弹毁伤数据的再分析

二.预制破片弹毁伤数据的再分析MC方法：在炸高2m-4m、落角45°-90°范围内，蒙特卡洛生产大量组合，代入回归模型得到预测毁伤面积后，统计超过2000平米的占比。importnumpy

asnp

data=np.array([

[-1,-1,1202.44],

[-1,-1,1509.92],

[-1,1,3171.77],

[-1,1,3486.43],

[1,-1,2742.72],

[1,-1,2423.20],

[1,1,2501.86],

[1,1,3155.33]])

x1_all=data[:,0]

x2_all=data[:,1]

y_all=data[:,2]

X=np.column_stack([np.ones(len(y_all),x1_all,x2_all,x1_all*x2_all])beta_hat=np.linalg.inv(X.T@X)@X.T@y_alln_samples=1000000

x1_mc=np.random.uniform(-1,1,n_samples)

x2_mc=np.random.uniform(-1,1,n_samples)beta_0,beta_1,beta_2,beta_12=beta_hat

y_pred=beta_0+beta_1*x1_mc+beta_2*x2_mc+beta_12*x1_mc*x2_mc

threshold=2000prob=np.mean(y_pred>threshold)

print(f"\n蒙特卡洛模拟次数:{n_samples}")print(f"预测y>{threshold}的概率:{prob:.4f}")约89.67%的条件下毁伤面积超过2000平方米！钻地弹侵彻深度数据的再