python数据分析— ch12-13 python参数估计与假设检验

python数据分析— ch12-13 python参数估计与假设检验

• 1. Ch12–python 参数估计
• • 1.1 参数估计与置信区间的含义及函数版
  • • 1.1.1 参数估计与置信区间的含义
    • 1.1.2 参数估计函数版
    • 1.1.3 参数估计函数版
  • 1.2 Python单正态总体均值区间估计
  • • 1.2.1 方差

      σ

      2

      \sigma^2

      σ2已知

    • 1.2.2 方差

      σ

      2

      \sigma^2

      σ2未知

  • 1.3 Python单正态总体方差区间估计
  • 1.4 Python双正态总体均值差区间估计
  • • 1.4.1 方差

      σ

      2

      \sigma^2

      σ2已知

    • 1.4.2 方差

      σ

      2

      \sigma^2

      σ2未知，但相等

      σ

      1

      2

      =

      σ

      2

      2

      \sigma_{1}^2=\sigma_{2}^2

      σ12​=σ22​

  • 1.5 两方差都未知时两均值差的置信区间
  • 1.6 Python双正态总体方差比区间估计
• 2. Ch13–python 假设检验
• • 2.1 Python单个样本t检验
  • 2.2 Python两个独立样本t检验
  • 2.3 Python配对样本t检验
  • 2.4 两独立样本 vs 配对样本
  • • 2.4.1 两个独立样本t检验
    • 2.4.2 配对样本t检验
    • 2.4.3 总结
  • 2.5 Python单样本方差假设检验
  • 2.6 Python双样本方差假设检验

1. Ch12–python 参数估计

1.1 参数估计与置信区间的含义及函数版

1.1.1 参数估计与置信区间的含义

参数估计和置信区间是统计学中的两个重要概念，它们在数据分析和推断统计中扮演着关键角色。下面分别解释这两个概念：

1.参数估计（Parameter Estimation）

参数估计是指使用样本数据来估计总体参数的过程。在统计学中，总体参数是指描述整个数据集的特征值，如总体均值（μ）、总体方差（σ²）等。由于我们通常无法获得整个总体的数据，因此需要从总体中抽取一部分样本，然后通过样本数据来推断总体参数。

参数估计有两种主要类型：

• 点估计（Point Estimation）：这是最基本的估计方法，它提供了一个单一的数值作为总体参数的估计。例如，样本均值（(\bar{x})）通常被用作总体均值（μ）的点估计。
  

  x

  ˉ

  =

  ∑

  i

  =

  1

  n

  x

  i

  n

  \bar{x} = \frac{\sum_{i=1}^{n} x_i}{n}

  xˉ=n∑i=1n​xi​​
  其中，

  x

  ˉ

  \bar{x}

  xˉ表示均值，

  x

  i

  x_i

  xi​表示第 𝑖个数值，而 𝑛是数值的总数,即样本数。

• 区间估计（Interval Estimation）：这种方法提供了一个数值范围，这个范围以一定的概率包含总体参数。区间估计通常比点估计提供更多信息，因为它包含了不确定性的度量。

例12-1：对某个篮球运动员记录其在某一次比赛中投篮命中与否，观测数据如下：

1 1 0 1 0 0 1 0 1 1 1 0 1 1 0 1

0 0 1 0 1 0 1 0 0 1 1 0 1 1 0 1

编写Python程序估计这个篮球运动员投篮的成败比。

import numpy as np
x=[1,1,0,1,0,0,1,0,1,1,1,0,1,1,0,1,0,0,1,0,1,0,1,0,0,1,1,0,1,1,0,1]
theta=np.mean(x)
h=theta/(1-theta)
print ('h=',round(h,4))

output
h= 1.2857

2.置信区间（Confidence Interval）

置信区间是区间估计的一种形式，它表示我们对总体参数估计的可信程度。一个置信区间通常表示为：估计值 ± 误差范围，其中误差范围被称为置信区间的“半宽”。数学表达式如下：

置信区间

=

x

^

±

z

α

/

2

×

S

E

(

x

^

)

ext{置信区间} = \hat{x} \pm z_{\alpha/2} imes SE(\hat{x})

置信区间=x^±zα/2​×SE(x^)

其中，$ z_{\alpha/2} $ 是标准正态分布的 $ \alpha/2 $ 分位数，$ SE(\hat{x})$ 是

x

^

\hat{x}

x^ 的标准误差。
置信区间的关键组成部分包括：

• 置信水平（Confidence Level）：这是一个百分比，表示我们有多大的信心认为置信区间包含了总体参数。常见的置信水平有90%、95%和99%。

• 置信区间的宽度：这取决于样本大小、样本变异性以及所选的置信水平。较大的样本大小和较低的变异性通常会导致更窄的置信区间。

• 总体参数：置信区间试图估计的总体特征，如均值、比例或差异。

例如，如果我们说“总体均值的95%置信区间是20到30”，这意味着我们有95%的信心认为总体均值位于20到30的范围内，基于我们从样本中得到的信息。

例12-2：假设一位投资分析师从股权基金中选取了一个随机样本，并计算出了平均的夏普比率。样本的容量为100，并且平均的夏普比率为0.45。该样本具有的标准差为0.300。利用一个基于标准正态分布的临界值，编写python程序计算并解释所有股权基金总体均值的90%置信区间。（

z

0.05

=

1.65

z_{0.05}=1.65

z0.05​=1.65）

import numpy as np
up=0.45+1.65*(0.300/np.sqrt(100))
low=0.45-1.65*(0.300/np.sqrt(100))
print(f"置信区间为：[{round(low,4)},{round(up,4)}]")

output
置信区间为：[0.4005,0.4995]

1.1.2 参数估计函数版

import scipy.stats as ss
import math
#confidence_interval
def norm1_ci_4mean(n, x_bar, alpha, sigma=None, s=None):
"""
计算单正态总体均值的置信区间。
参数:
n : int
样本大小。
x_bar : float
样本均值。
alpha : float
置信水平，例如0.05表示95%置信区间。
sigma : float, optional
总体标准差，如果为None，则假定总体方差未知。
s : float, optional
样本标准差，如果为None且sigma也为None，则假定总体方差未知。
返回:
tuple
均值的置信区间下限和上限。
"""
# 计算置信区间的临界值
if sigma is None and s is not None:
if n<30:
# 总体方差未知，使用t分布
critical = ss.t.ppf(1 - alpha / 2, df=n - 1)
else:
critical = ss.norm.ppf(1 - alpha / 2)
elif sigma is not None and s is None:
# 总体方差已知，使用正态分布
critical = ss.norm.ppf(1 - alpha / 2)
else:
raise ValueError("参数sigma和s不能同时为None或非None。")
print(critical)
# 计算置信区间
if sigma is None:
# 使用样本标准差和t分布
margin_of_error = critical * (s / math.sqrt(n))
else:
# 使用总体标准差和正态分布
margin_of_error = critical * (sigma / math.sqrt(n))
low_bound = x_bar - margin_of_error
up_bound = x_bar + margin_of_error
return [round(low_bound, 4), round(up_bound, 4)]
def norm1_ci_4var(n, var, alpha):
"""
计算单正态总体方差的置信区间。
参数:
n : int
样本大小。
var : float
样本方差。
alpha : float
置信水平，例如0.05表示95%置信区间。
返回:
list
方差的置信区间下限和上限。
"""
# 计算置信区间的临界值
low_critical = ss.chi2.ppf(1 - alpha/2, df=n - 1)
up_critical = ss.chi2.ppf(alpha/2, df=n - 1)
low_bound = (n-1)*var/low_critical
up_bound = (n-1)*var/up_critical
return [round(low_bound, 4), round(up_bound, 4)]
def norm2_ci_4difference_mean(x1_bar,x2_bar,n1,n2,alpha,s1_2=None,s2_2=None,sigma1_2=None,sigma2_2=None):
"""
计算两个独立正态总体均值之差的置信区间。
参数:
mu1 : float
第一个样本的均值。
s1 : float
第一个样本的标准差。
n1 : int
第一个样本的大小。
mu1 : float
第二个样本的均值。
s2 : float
第二个样本的标准差。
n2 : int
第二个样本的大小。
alpha : float
置信水平，例如0.05表示95%置信区间。
equal_var : bool, optional
如果为True，则假定两个总体方差相等，使用z统计量；如果为False，则假定方差不相等，使用t统计量。
返回:
list
均值差值的置信区间下限和上限。
"""
critical_value = 0
std_error = 0
# 计算均值差值的标准误差
if sigma1_2 is not None and sigma2_2 is not None:#两总体方差已知
std_error = math.sqrt(sigma1_2 / n1 + sigma2_2 / n2)
critical_value = ss.norm.ppf(1 - alpha/2)
if s1_2 is not None and s2_2 is not None:#两总体方差未知但相等
s_2 = ((n1-1)*s1_2 + (n2-1)*s2_2)/((n1-1)+(n2-1))
std_error = math.sqrt((1/n1 + 1/n2)*s_2)
critical_value = ss.t.ppf(1 - alpha/2,df = n1+n2-2)
# 计算置信区间
margin_of_error = critical_value * std_error
low_bound = x1_bar-x2_bar - margin_of_error
up_bound = x1_bar-x2_bar + margin_of_error
return [round(low_bound, 4), round(up_bound, 4)]
def norm2_ci_4var(s1_2,s2_2,n1,n2,alpha):
"""
计算两个独立正态总体均值之差的置信区间。
参数:
s1 : float
第一个样本的标准差。
n1 : int
第一个样本的大小。
s2 : float
第二个样本的标准差。
n2 : int
第二个样本的大小。
alpha : float
置信水平，例如0.05表示95%置信区间。
返回:
list
均值差值的置信区间下限和上限。
"""
low_critical_value = 1/ss.f.ppf(1-alpha/2, n1-1, n2-1)
up_critical_value = 1/ss.f.ppf(alpha/2, n1-1, n2-1)
# 计算均值差值的标准误差
low_bound = s1_2/s2_2*low_critical_value
up_bound = s1_2/s2_2*up_critical_value
return [round(low_bound, 4), round(up_bound, 4)]

# 示例使用
n = 100
x_bar = 0.45
alpha = 0.1
sigma = None  # 假设总体方差未知
s = 0.3      # 样本标准差
# 调用函数
ci = norm1_confidence_interval_4mean(n, x_bar, alpha, sigma, s)
print(f"{(1-alpha)*100}%置信区间为：{ci}")

output
90.0%置信区间为：[0.4007, 0.4993]

1.1.3 参数估计函数版

from scipy import stats
import math
from typing import Tuple, Optional
class ParameterEstimation:
@staticmethod
def norm1_ci_for_mean(n: int, x_bar: float, alpha: float, sigma: Optional[float] = None, s: Optional[float] = None) -> Tuple[float, float]:
"""
计算单正态总体均值的置信区间。
"""
if sigma is not None and s is not None:
raise ValueError("Provide either sigma (known population standard deviation) or s (sample standard deviation), not both.")
critical_value = stats.t.ppf(1 - alpha / 2, df=n - 1) if s is not None else stats.norm.ppf(1 - alpha / 2)
margin_of_error = critical_value * (s / math.sqrt(n) if s is not None else sigma)
return (round(x_bar - margin_of_error, 4), round(x_bar + margin_of_error, 4))
@staticmethod
def norm1_ci_for_var(n: int, sample_var: float, alpha: float) -> Tuple[float, float]:
"""
计算单正态总体方差的置信区间。
"""
critical_value_low = stats.chi2.ppf(1 - alpha / 2, df=n - 1)
critical_value_high = stats.chi2.ppf(alpha / 2, df=n - 1)
lower_bound = (n - 1) * sample_var / critical_value_low
upper_bound = (n - 1) * sample_var / critical_value_high
return (round(lower_bound, 4), round(upper_bound, 4))
@staticmethod
def norm2_ci_for_difference_of_means(x_bar1: float, x_bar2: float, n1: int, n2: int, alpha: float, s1: Optional[float] = None, s2: Optional[float] = None, sigma1: Optional[float] = None, sigma2: Optional[float] = None) -> Tuple[float, float]:
"""
计算两个独立正态总体均值之差的置信区间。
"""
if (sigma1 is not None and sigma2 is not None) and (s1 is not None and s2 is not None):
raise ValueError("For known population variances, do not provide sample standard deviations.")
std_error = math.sqrt(sigma1 / n1 + sigma2 / n2) if sigma1 is not None and sigma2 is not None else math.sqrt(((n1 - 1) * s1 + (n2 - 1) * s2) / (n1 + n2 - 2))
critical_value = stats.t.ppf(1 - alpha / 2, df=n1 + n2 - 2) if s1 is not None and s2 is not None else stats.norm.ppf(1 - alpha / 2)
margin_of_error = critical_value * std_error
difference = x_bar1 - x_bar2
return (round(difference - margin_of_error, 4), round(difference + margin_of_error, 4))
@staticmethod
def norm2_ci_for_var_ratio(s1: float, s2: float, n1: int, n2: int, alpha: float) -> Tuple[float, float]:
"""
计算两个独立正态总体方差比的置信区间。
"""
critical_value_low = 1 / stats.f.ppf(1 - alpha / 2, dfn=n1 - 1, dfd=n2 - 1)
critical_value_high = 1 / stats.f.ppf(alpha / 2, dfn=n1 - 1, dfd=n2 - 1)
lower_bound = s1 / s2 * critical_value_low
upper_bound = s1 / s2 * critical_value_high
return (round(lower_bound, 4), round(upper_bound, 4))

# 示例使用
try:
ci_mean = ParameterEstimation.norm1_ci_for_mean(n=100, x_bar=50, alpha=0.05, s=5)
print(f"Confidence interval for mean: {ci_mean}")
except ValueError as e:
print(e)

output
Confidence interval for mean: (49.0079, 50.9921)

• 总结

参数估计是使用样本数据来推断总体特征的过程，而置信区间是参数估计的一种形式，它提供了一个数值范围，表示我们对总体参数估计的可信程度。置信区间的宽度和所选的置信水平反映了估计的不确定性。在实际应用中，选择合适的置信水平和理解置信区间的含义对于做出准确的统计推断至关重要。

1.2 Python单正态总体均值区间估计

1.2.1 方差

σ

2

\sigma^2

σ2已知

置信区间

=

x

ˉ

±

Z

α

2

σ

n

=

(

x

ˉ

−

Z

α

2

σ

n

,

x

ˉ

+

Z

1

−

α

2

σ

n

)

ext{置信区间} = \bar{x} \pm Z_{\frac{\alpha}{2}} \frac{\sigma}{\sqrt{n}}=(\bar{x} – Z_{\frac{\alpha}{2}} \frac{\sigma}{\sqrt{n}},\bar{x} + Z_{1-\frac{\alpha}{2}} \frac{\sigma}{\sqrt{n}})

置信区间=xˉ±Z2α​​n

​σ​=(xˉ−Z2α​​n

​σ​,xˉ+Z1−2α​​n

​σ​)

例12-3：某车间生产的滚珠直径X服从正态分布N(u,0.6)。现从某天的产品中抽取6个，测得直径如下(单位：mm)：

14.6，15.1，14.9，14.8，15.2，15.1

编写python程序试求置信度为95%的平均直径置信区间。

import numpy as np
import scipy.stats as ss
n = 6; p = 0.025; sigma = np.sqrt(0.6)
x=[14.6,15.1,14.9,14.8,15.2,15.1]
xbar=np.mean(x)
low = xbar - ss.norm.ppf(q = 1 - p) * (sigma / np.sqrt(n))
up = xbar + ss.norm.ppf(q = 1 - p) * (sigma / np.sqrt(n))
print ('low=',round(low,4))
print ('up=',round(up,4))

output
low= 14.3302
up= 15.5698

import numpy as np
n = 6
x=[14.6,15.1,14.9,14.8,15.2,15.1]
x_bar = np.mean(x)
alpha = 0.05
sigma = np.sqrt(0.6)  # 总体标准差
s = None      # 样本标准差
# 调用函数
ci = norm1_confidence_interval_4mean(n, x_bar, alpha, sigma, s)
print(f"{(1-alpha)*100}%置信区间为：{ci}")

output
1.959963984540054
95.0%置信区间为：[14.3302, 15.5698]

1.2.2 方差

σ

2

\sigma^2

σ2未知

置信区间

=

x

ˉ

±

t

α

2

(

n

−

1

)

S

n

=

(

x

ˉ

−

t

α

2

(

n

−

1

)

S

n

,

x

ˉ

+

t

1

−

α

2

(

n

−

1

)

S

n

)

ext{置信区间} = \bar{x} \pm t_{\frac{\alpha}{2}}(n-1) \frac{S}{\sqrt{n}}=(\bar{x} – t_{\frac{\alpha}{2}}(n-1) \frac{S}{\sqrt{n}} , \bar{x} + t_{1-\frac{\alpha}{2}}(n-1) \frac{S}{\sqrt{n}})

置信区间=xˉ±t2α​​(n−1)n

​S​=(xˉ−t2α​​(n−1)n

​S​,xˉ+t1−2α​​(n−1)n

​S​)

例12-4：某糖厂自动包装机装糖，设各包重量服从正态分布N(u,

σ

2

\sigma^2

σ2)。 某日开工后测得9包重量为(单位：kg)：

99.10，98.7，100.5，101.2，98.10，99.7，99.5，102.1，100.5

编写python程序试求置信度为95%的平均直径置信区间。

import numpy as np
import scipy.stats as ss
from scipy.stats import t
n = 9; p = 0.025; s = np.sqrt(1.47)
x=[99.3,98.7,100.5,101.2,98.3,99.7,99.5,102.1,100.5]
xbar=np.mean(x)
low = xbar - ss.t.ppf(1-p,n-1) * (s / np.sqrt(n))
up = xbar + ss.t.ppf(1-p,n-1) * (s / np.sqrt(n))
print ('low=',round(low,4))
print ('up=',round(up,4))

output
low= 99.0458
up= 100.9097

import numpy as np
x=[99.3,98.7,100.5,101.2,98.3,99.7,99.5,102.1,100.5]
# 调用函数
alpha=0.05
ci = norm1_confidence_interval_4mean(n=9, x_bar=np.mean(x), alpha=0.05, sigma=None, s=np.sqrt(1.47))
print(f"{(1-alpha)*100}%置信区间为：{ci}")

output
2.3060041350333704
95.0%置信区间为：[99.0458, 100.9097]

from scipy import stats
import numpy as np
x=[99.3,98.7,100.5,101.2,98.3,99.7,99.5,102.1,100.5]
stats.t.interval(0.95,len(x)-1,np.mean(x),stats.sem(x))

output
(99.04599342616191, 100.90956212939363)

1.3 Python单正态总体方差区间估计

此时虽然也可以就均值是否已知分两种情况讨论方差的区间估计，但在实际中u已知的情形是极为罕见的，所以只讨论在u未知的情况下，方差

σ

2

\sigma^2

σ2的置信区间。

例12-5 从某车间加工的同类零件中抽取了16件，测得零件的平均长度为12.8厘米，方差为0.0023。假设零件的长度服从正态分布，试求总体方差及标准差的置信区间(置信度为95%)。

from scipy.stats import chi2
n=16;sq=0.0023;p=0.025
low = ((n-1)*sq)/ chi2.ppf(1-p, n-1)
up = ((n-1)*sq)/ chi2.ppf(p, n-1)
print ('low=',round(low,4))
print ('up=',round(up,4))

output
low= 0.0013
up= 0.0055

norm1_confidence_interval_4var(n=16, var=0.0023, alpha=0.05)

output
[0.0013, 0.0055]

1.4 Python双正态总体均值差区间估计

1.4.1 方差

σ

2

\sigma^2

σ2已知

假设

σ

1

2

,

σ

2

2

\sigma_1^2,\sigma_2^2

σ12​,σ22​都已知，求

μ

1

−

μ

2

\mu_1-\mu_2

μ1​−μ2​的置信度为

1

−

α

1 – \alpha

1−α的置信区间

由于

X

ˉ

∼

N

(

μ

1

,

σ

1

2

)

\bar{X} \sim N(\mu_1, \sigma_1^2)

Xˉ∼N(μ1​,σ12​),

Y

ˉ

∼

N

(

μ

2

,

σ

2

2

)

\bar{Y} \sim N(\mu_2, \sigma_2^2)

Yˉ∼N(μ2​,σ22​)，并且两者独立，得到：

X

ˉ

−

Y

ˉ

∼

N

(

μ

1

−

μ

2

,

σ

1

2

n

1

+

σ

2

2

n

2

)

\bar{X} – \bar{Y} \sim N(\mu_1 – \mu_2, \frac{\sigma_1^2}{n_1} + \frac{\sigma_2^2}{n_2})

Xˉ−Yˉ∼N(μ1​−μ2​,n1​σ12​​+n2​σ22​​)

因此有：

Z

=

(

X

ˉ

−

Y

ˉ

)

−

(

μ

1

−

μ

2

)

σ

1

2

n

1

+

σ

2

2

n

2

∼

N

(

0

,

1

)

Z=\frac{(\bar{X} – \bar{Y}) – (\mu_1 – \mu_2)}{\sqrt{\frac{\sigma_1^2}{n_1} + \frac{\sigma_2^2}{n_2}}}\sim N(0,1)

Z=n1​σ12​​+n2​σ22​​

​(Xˉ−Yˉ)−(μ1​−μ2​)​∼N(0,1)

由

P

(

−

z

α

/

2

<

Z

<

z

1

−

α

/

2

)

=

1

−

α

P(-z_{\alpha/2} < Z < z_{1-\alpha/2}) = 1 – \alpha

P(−zα/2​<Z<z1−α/2​)=1−α 即得：

P

(

X

ˉ

−

Y

ˉ

−

z

α

/

2

σ

1

2

n

1

+

σ

2

2

n

2

<

μ

1

−

μ

2

<

X

ˉ

−

Y

ˉ

+

z

α

/

2

σ

1

2

n

1

+

σ

2

2

n

2

)

=

1

−

α

P\left(\bar{X} – \bar{Y} – z_{\alpha/2}\sqrt{\frac{\sigma_1^2}{n_1} + \frac{\sigma_2^2}{n_2}}< \mu_1 – \mu_2 < \bar{X} – \bar{Y} + z_{\alpha/2}\sqrt{\frac{\sigma_1^2}{n_1} + \frac{\sigma_2^2}{n_2}}\right) = 1 – \alpha

P

​Xˉ−Yˉ−zα/2​n1​σ12​​+n2​σ22​​

​<μ1​−μ2​<Xˉ−Yˉ+zα/2​n1​σ12​​+n2​σ22​​

​

​=1−α

所以两均值差的置信区间为：

(

X

ˉ

−

Y

ˉ

−

z

α

/

2

σ

1

2

n

1

+

σ

2

2

n

2

,

X

ˉ

−

Y

ˉ

+

z

α

/

2

σ

1

2

n

1

+

σ

2

2

n

2

)

\left(\bar{X} – \bar{Y} – z_{\alpha/2}\sqrt{\frac{\sigma_1^2}{n_1} + \frac{\sigma_2^2}{n_2}}, \bar{X} – \bar{Y} + z_{\alpha/2}\sqrt{\frac{\sigma_1^2}{n_1} + \frac{\sigma_2^2}{n_2}}\right)

​Xˉ−Yˉ−zα/2​n1​σ12​​+n2​σ22​​

​,Xˉ−Yˉ+zα/2​n1​σ12​​+n2​σ22​​

​

​

同理可求得，两均值差的置信度为 $ 1 – \alpha $ 的置信上限为：

X

ˉ

−

Y

ˉ

+

Z

1

−

α

/

2

σ

1

2

n

1

+

σ

2

2

n

2

\bar{X} – \bar{Y} + Z_{1-\alpha/2}\sqrt{\frac{\sigma_1^2}{n_1} + \frac{\sigma_2^2}{n_2}}

Xˉ−Yˉ+Z1−α/2​n1​σ12​​+n2​σ22​​

​

两均值差的置信度为 $ 1 – \alpha $ 的置信下限为：

X

ˉ

−

Y

ˉ

−

z

1

−

α

/

2

σ

1

2

n

1

+

σ

2

2

n

2

\bar{X} – \bar{Y} – z_{1-\alpha/2}\sqrt{\frac{\sigma_1^2}{n_1} + \frac{\sigma_2^2}{n_2}}

Xˉ−Yˉ−z1−α/2​n1​σ12​​+n2​σ22​​

​

例12-6：为比较两种农产品的产量，选择18块条件相似的试验田，采用相同的耕作方法做实验，结果播种甲品种的8块试验田的单位面积产量和播种乙品种的10块试验田的单
位面积产量分别如表12-2所示。

两种农产品的产量如下：

假定每个品种的单位面积产量均服从正态分布，甲品种产量的方差为2140，乙品种产量的方差为3250，试求这两个品种平均面积产量差的置信区间（取a=0.05)。

import numpy as np
import scipy.stats as ss
x=[628,583,510,554,612,523,530,615]
y=[535,433,398,470,567,480,498,560,503,426]
n1=len(x)
n2=len(y)
x_bar=np.mean(x)
y_bar=np.mean(y)
sigma1_2=2140
sigma2_2=3250
alpha = 0.05
low = x_bar - y_bar-ss.norm.ppf(q = 1 - alpha/2) * np.sqrt(sigma1_2/n1+sigma2_2/n2)
up = x_bar - y_bar+ss.norm.ppf(q = 1 - alpha/2) * np.sqrt(sigma1_2/n1+sigma2_2/n2)
print ('low=',round(low,4))
print ('up=',round(up,4))

output
low= 34.6669
up= 130.0831

ci = norm2_ci_4difference_mean(x_bar,y_bar,n1,n2,alpha,s1_2=None,s2_2=None,sigma1_2=sigma1_2,sigma2_2=sigma2_2)
print(f"{(1-alpha)*100}%置信区间为：{ci}")

output
95.0%置信区间为：[34.6669, 130.0831]

1.4.2 方差

σ

2

\sigma^2

σ2未知，但相等

σ

1

2

=

σ

2

2

\sigma_{1}^2=\sigma_{2}^2

σ12​=σ22​

1.5 两方差都未知时两均值差的置信区间

设两方差均未知，但

σ

1

2

=

σ

2

2

\sigma_1^2 = \sigma_2^2

σ12​=σ22​。此时由于

Z

=

(

X

ˉ

−

Y

ˉ

)

−

(

μ

1

−

μ

2

)

σ

1

2

n

1

+

σ

2

2

n

2

∼

N

(

0

,

1

)

Z=\frac{(\bar{X} – \bar{Y}) – (\mu_1 – \mu_2)}{\sqrt{\frac{\sigma_1^2}{n_1} + \frac{\sigma_2^2}{n_2}}}\sim N(0,1)

Z=n1​σ12​​+n2​σ22​​

​(Xˉ−Yˉ)−(μ1​−μ2​)​∼N(0,1) ,

(

n

1

−

1

)

S

1

2

σ

2

∼

χ

2

(

n

1

−

1

)

\frac{(n_1-1)S1^2}{\sigma^2} \sim \chi^2(n_1-1)

σ2(n1​−1)S12​∼χ2(n1​−1),

(

n

2

−

1

)

S

2

2

σ

2

∼

χ

2

(

n

2

−

1

)

\frac{(n_2-1)S2^2}{\sigma^2} \sim \chi^2(n_2-1)

σ2(n2​−1)S22​∼χ2(n2​−1)

所以

(

n

1

−

1

)

S

1

2

σ

1

2

+

(

n

2

−

1

)

S

2

2

σ

2

2

∼

χ

2

(

n

1

+

n

2

−

2

)

\frac{(n_1-1)S_1^2}{\sigma_1^2} + \frac{(n_2-1)S_2^2}{\sigma_2^2} \sim \chi^2(n_1+n_2-2)

σ12​(n1​−1)S12​​+σ22​(n2​−1)S22​​∼χ2(n1​+n2​−2)。

由此可得：

T

=

(

X

ˉ

−

Y

ˉ

−

(

μ

1

−

μ

2

)

(

1

n

1

+

1

n

2

)

S

2

∼

t

(

n

1

+

n

2

−

1

)

T=\frac{(\bar{X} – \bar{Y} – (\mu_1 – \mu_2)}{\sqrt{(\frac{1}{n_1}+\frac{1}{n_2})S^2}} \sim t(n_1+n_2-1)

T=(n1​1​+n2​1​)S2

​(Xˉ−Yˉ−(μ1​−μ2​)​∼t(n1​+n2​−1)。

其中，

S

2

=

(

n

1

−

1

)

S

1

2

+

(

n

2

−

1

)

S

2

2

(

n

1

−

1

)

+

(

n

2

−

1

)

S^2 = \frac{(n_1-1)S_1^2+(n_2-1)S_2^2}{(n_1-1)+(n_2-1)}

S2=(n1​−1)+(n2​−1)(n1​−1)S12​+(n2​−1)S22​​,

S

1

S_1

S1​和

S

2

S_2

S2​分别是样本1和样本2的标准差。

同样由

P

(

−

t

1

−

α

/

2

(

n

1

+

n

2

−

2

)

<

T

<

t

1

−

α

/

2

(

n

1

+

n

2

−

2

)

)

=

1

−

α

P\left(-t_{1-\alpha/2}(n_1+n_2-2) < T < t_{1-\alpha/2}(n_1+n_2-2)\right) = 1 – \alpha

P(−t1−α/2​(n1​+n2​−2)<T<t1−α/2​(n1​+n2​−2))=1−α，

解不等式即得两均值差的置信度为

1

−

α

1 – \alpha

1−α 的置信区间为：

(

X

ˉ

−

Y

ˉ

±

t

1

−

α

/

2

(

n

1

+

n

2

−

2

)

⋅

(

1

n

1

+

1

n

2

)

S

2

)

\left(\bar{X} – \bar{Y} \pm t_{1-\alpha/2}(n_1+n_2-2) \cdot \sqrt{(\frac{1}{n_1}+\frac{1}{n_2})S^2}\right)

(Xˉ−Yˉ±t1−α/2​(n1​+n2​−2)⋅(n1​1​+n2​1​)S2

​)

同理可求得两均值差的置信度为

1

−

α

1 – \alpha

1−α的置信上限为：

(

X

ˉ

−

Y

ˉ

+

t

1

−

α

/

2

(

n

1

+

n

2

−

2

)

⋅

(

1

n

1

+

1

n

2

)

S

2

)

\left(\bar{X} – \bar{Y} + t_{1-\alpha/2}(n_1+n_2-2) \cdot \sqrt{(\frac{1}{n_1}+\frac{1}{n_2})S^2}\right)

(Xˉ−Yˉ+t1−α/2​(n1​+n2​−2)⋅(n1​1​+n2​1​)S2

​)

两均值差的置信度为

1

−

α

1 – \alpha

1−α 的置信下限为：

(

X

ˉ

−

Y

ˉ

−

t

1

−

α

/

2

(

n

1

+

n

2

−

2

)

⋅

(

1

n

1

+

1

n

2

)

S

2

)

\left(\bar{X} – \bar{Y} – t_{1-\alpha/2}(n_1+n_2-2) \cdot \sqrt{(\frac{1}{n_1}+\frac{1}{n_2})S^2}\right)

(Xˉ−Yˉ−t1−α/2​(n1​+n2​−2)⋅(n1​1​+n2​1​)S2

​)

在上一例题12-6中，如果不知道两种品种产量的方差，但已知两者相同，试求这两个品种平均面积产量差的置信区间（取a=0.05)。

import numpy as np
import scipy.stats as ss
x=[628,583,510,554,612,523,530,615]
y=[535,433,398,470,567,480,498,560,503,426]
n1=len(x)
n2=len(y)
s1_2=np.var(x)
s2_2=np.var(y)
x_bar=np.mean(x)
y_bar=np.mean(y)
alpha = 0.05
sq=((n1-1)*s1_2+(n2-1)*s2_2)/(n1-1+n2-1)
low = xbar - ybar-ss.t.ppf(1-alpha/2,n1+n2-2)*np.sqrt(sq*(1/n1+1/n2))
up = xbar - ybar+ss.t.ppf(1-alpha/2,n1+n2-2)*np.sqrt(sq*(1/n1+1/n2))
print ('low=',round(low,4))
print ('up=',round(up,4))

output
low= 32.4209
up= 132.3291

ci = norm2_ci_4difference_mean(x_bar,y_bar,n1,n2,alpha,s1_2,s2_2,sigma1_2=None,sigma2_2=None)
print(f"{(1-alpha)*100}%置信区间为：{ci}")

output
95.0%置信区间为：[32.4209, 132.3291]

1.6 Python双正态总体方差比区间估计

此时虽然可以就均值是否已知分两种情况讨论方差的区间估计，但在实际中，

μ

\mu

μ已知的情形是极为罕见的，所以只讨论在没有位置的情况下，方差

σ

2

\sigma^2

σ2的置信区间。

由于

(

n

1

−

1

)

S

1

2

σ

2

∼

χ

2

(

n

1

−

1

)

\frac{(n_1-1)S_1^2}{\sigma^2} \sim \chi^2(n_1-1)

σ2(n1​−1)S12​​∼χ2(n1​−1),