SCAU统计学实验

实验一

描述：

第一次实验有 $n$ 个人抛掷硬币，若有一人的硬币与其它 $n-1$ 人的不同，则游戏结束。

问：第一次就结束的概率是多少？求游戏结束时抛掷硬币的平均次数。请结合概率的统计定义，编程仿真实现并求解以上问题。

已知： $n=6$ 时，概率为 $0.1875$ （仅用于仿真比较）

分析:

$$p(\text{第一次结束概率})=2\times C_{n}^{1}\cdot0.5^{\prime}\cdot(1-0.5)^{n-1}$$

抛掷硬币的平均次数，即求抛掷硬币次数的数学期望，当 $n$ 足够大时，期望 $\approx$ 平均次数。

$$E(x)=\sum_{i=1}^{n}i\cdot p(\text{第 i-1 轮未结束})\cdot p(\text{第 i 轮结束})$$

由于每次抛掷硬币都是独立的，所以 $p(\text{第 i 轮结束}) = p(\text{第一次结束概率})$。

以上思路是基于概率学的推导计算，也可以采用蒙特卡罗法，用 random 模拟每次的结果，从而计算出结果。

实现代码:

from math import *


def p_end(n):
    return 2 * comb(n, 1) * pow(0.5, n)


def train(n, epoch):
    p_i_exit = p_end(n)  # 第i轮退出的概率
    p_i_1_run = 1  # 前i-1轮未退出的概率
    Ex = 1  # 期望
    for i in range(1, epoch):
        p_i_1_run *= (1 - p_i_exit)
        Ex += i * p_i_exit * p_i_1_run
    return Ex


if __name__ == "__main__":
    n = 6
    print("第一次结束概率：", p_end(n))
    print("抛掷硬币的平均次数：{:.6f}".format(train(n, 1000000)))

输出结果：

第一次结束概率： 0.1875
抛掷硬币的平均次数：5.333333

实验二

描述：验证中心极限定理。

分析:

step 1：首先生成一个总体（约10000个数据），计算出总体的均值与方差

step 2：指定一个自然数 $n$（$n$可分别取2，3，5，15，25，30 …）

step 3：从总体中抽取 $n$ 个数据组成一个样本，计算出样本的均值

step 4：重复step 3 $m$ 次（例如：$m = 1000$）

step 5：计算 $m$ 个样本均值的平均值和方差

step 6：验证 $m$ 个样本均值的平均值和方差与总体的均值、方差的数理关系

step 7：重复step 2到step 6的过程， $n$ 取不同的对应值

step 8：生成不一样的总体数据，重复step 1到step 7的过程

实现代码:

import numpy as np


def test2():
    # step 1：生成数据，求均值、方差
    size = 10000
    data = np.random.normal(loc=100, scale=10, size=size)  # 生成均值为100，标准差为10的总体数据
    tot_mean = np.mean(data)  # 总体均值
    tot_var = np.var(data)  # 总体方差

    # step 2：指定批次大小、抽取次数
    batch_size = [2, 3, 5, 15, 25, 30, 50]
    epoch = 1000

    # step 3：抽取 batch_size 个数据，计算出样本的均值 
    # step 4：重复step 3 epoch 次

    for n in batch_size:
        print("批次大小为", n, "时：")
        sample_means = []  # 样本均值
        for i in range(epoch):
            sample = np.random.choice(data, size=n, replace=True)
            sample_means.append(np.mean(sample))
        # step 5：计算 epoch 个样本均值的平均值和方差
        sample_means_mean = np.mean(sample_means)
        sample_means_var = np.var(sample_means)
        # step 6：验证 epoch 个 样本均值的平均值和方差 与 总体的均值、方差 的关系
        print(f"    总体均值: {tot_mean:.2f}, 样本均值的平均值: {sample_means_mean:.2f}")
        print(f"    总体方差: {tot_var:.2f}, 样本均值的方差: {sample_means_var:.2f}")
        print()


if __name__ == "__main__":
    for _ in range(100):
        print(f"************第{_ + 1}次实验************")
        test2()

输出结果：

第1次实验
批次大小为 2 时：
总体均值: 100.09, 样本均值的平均值: 100.12
总体方差: 98.84, 样本均值的方差: 51.28

批次大小为 3 时：
总体均值: 100.09, 样本均值的平均值: 100.17
总体方差: 98.84, 样本均值的方差: 33.31

批次大小为 5 时：
总体均值: 100.09, 样本均值的平均值: 99.95
总体方差: 98.84, 样本均值的方差: 19.90

批次大小为 15 时：
总体均值: 100.09, 样本均值的平均值: 100.09
总体方差: 98.84, 样本均值的方差: 6.29

批次大小为 25 时：
总体均值: 100.09, 样本均值的平均值: 100.22
总体方差: 98.84, 样本均值的方差: 4.03

批次大小为 30 时：
总体均值: 100.09, 样本均值的平均值: 100.08
总体方差: 98.84, 样本均值的方差: 3.20

批次大小为 50 时：
总体均值: 100.09, 样本均值的平均值: 100.13
总体方差: 98.84, 样本均值的方差: 1.90

…

第100次实验
批次大小为 2 时：
总体均值: 100.02, 样本均值的平均值: 99.95
总体方差: 97.85, 样本均值的方差: 50.26

批次大小为 3 时：
总体均值: 100.02, 样本均值的平均值: 100.21
总体方差: 97.85, 样本均值的方差: 33.85

批次大小为 5 时：
总体均值: 100.02, 样本均值的平均值: 99.96
总体方差: 97.85, 样本均值的方差: 19.24

批次大小为 15 时：
总体均值: 100.02, 样本均值的平均值: 100.05
总体方差: 97.85, 样本均值的方差: 6.40

批次大小为 25 时：
总体均值: 100.02, 样本均值的平均值: 100.01
总体方差: 97.85, 样本均值的方差: 3.79

批次大小为 30 时：
总体均值: 100.02, 样本均值的平均值: 100.13
总体方差: 97.85, 样本均值的方差: 3.26

批次大小为 50 时：
总体均值: 100.02, 样本均值的平均值: 99.95
总体方差: 97.85, 样本均值的方差: 2.05

实验三

描述：

根据自己手机微信每月支付的数据计算出微信支付的月支付95%的置信区间
并画出近三年来月支付金额的分布图（为避免隐私问题，数据可以自己构造）

分析：

首先我们需要自己构造数据，并且近三年内月支付只有 3 * 12 = 36 个数据，属于小样本量，并且大概率不属于正态分布。

于是我们假设数据服从T分布，来计算置信区间。

这种方法适用于样本量较小或者数据不满足正态分布的情况。

实现代码:

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

# 解决matplotlib中文乱码
plt.rcParams['font.sans-serif'] = ['SimSun']

data = np.array([
    1326, 1032, 1445, 1612, 1455, 2176, 2172, 1508, 1593, 2356, 1713, 1871,
    1934, 1421, 2149, 1791, 1648, 1477, 1716, 1279, 1755, 1657, 1712, 1950,
    1577, 2005, 1787, 1768, 1837, 1823, 1325, 2145, 2050, 1901, 1908, 1744
])

# 计算样本均值、标准差、标准误差
mean = np.mean(data)
std = np.std(data)
se = std / np.sqrt(len(data))

# 计算置信区间
confidence = 0.95  # 置信度为95%
critical_value = stats.t.ppf((1 + confidence) / 2, len(data) - 1) * se  # 计算临界值
lower_limit = mean - critical_value  # 置信区间的下限
upper_limit = mean + critical_value  # 置信区间的上限

# 输出
print(f"微信月支付\n"
      f"置信区间：[{lower_limit:.4f}, {upper_limit:.4f}]\n"
      f"置信度：95%")

# 数据排序
data_sorted = np.sort(data)

# 统计频数
bin_size = 200
hist, bins = np.histogram(data_sorted, bins=np.arange(data_sorted.min(), data_sorted.max() + bin_size, bin_size))

# 绘制条形图
plt.bar(bins[:-1], hist, width=bin_size, align='edge', edgecolor='black', color='skyblue', label='频数')

# 绘制置信区间上下限
plt.axvline(lower_limit, color='red', linestyle='--', label='置信区间下限')
plt.axvline(upper_limit, color='green', linestyle='--', label='置信区间上限')

# 添加标题和标签
plt.title('月支付金额分布图')
plt.xlabel('金额')
plt.ylabel('频数')
plt.legend()

# 显示图形
plt.show()

输出结果：

微信月支付
置信区间：[1642.3460, 1836.4318]
置信度：95%

实验四

编码实现单个总体的假设检验

描述：

某地区小麦的一般生产水平为亩产 250 千克，其标准差为 30 千克。现用一种化肥进行试验， 从 25 个地块抽样，平均亩产量为 270 千克。

这种化肥是否使小麦明显增产（$\alpha=0.05$）?

分析：

总体平均值（$μ$）：250 kg

总体标准差（$\sigma^{2}$）：30kg

样本平均值（$\bar{x}$）：270 kg

样本大小（$n$）：25

我们已经知道小麦亩产的总体方差（$\sigma^{2}$），所以应该使用 $Z$ 检验。

需要检验使用化肥后平均亩产是否增产（$\mu>250$），即右侧检验。

$H_{0}:\mu\leq250$ $H_{1}:\mu>250$

$Z$ 检验统计量：

$$Z=\frac{\overline{x}-\mu}{\sigma/\sqrt{n}}$$

实现代码:

import math
import numpy as np
from scipy.stats import norm


def z_test(tot_mean, tot_var, sample_mean, n):
    # 计算标准误差和 Z 统计量
    SE = tot_var / math.sqrt(n)
    z_score = (sample_mean - tot_mean) / SE

    # 计算 p 值
    p_value = 1 - norm.cdf(np.abs(z_score))

    return z_score, p_value


if __name__ == "__main__":
    # 总体和样本参数
    tot_mean = 250  # 总体均值
    tot_var = 30  # 总体方差
    sample_mean = 270  # 样本均值
    n = 25  # 样本量
    alpha = 0.05

    z_score, p_value = z_test(tot_mean, tot_var, sample_mean, n)

    print(f"Z值: {z_score:.4f}")
    print(f"P值: {p_value:.4f}")

    if p_value < alpha:
        print("拒绝原假设，使用化肥后平均亩产增产")
    else:
        print("接受原假设，使用化肥后平均亩产未增产")

输出结果：

Z值: 3.3333
P值: 0.0004
拒绝原假设，使用化肥后平均亩产增产

经过$Z$ 检验后，拒绝原假设，该化肥能使小麦明显增产。

描述：

某种电子元件的寿命 $x$ 服从正态分布。现测得 16 只元件的寿命（单位：小时）如下：

159	280	101	212	224	379	179	264
222	362	168	250	149	260	485	170

是否有理由认为元件的平均寿命显著地大于225 小时（$\alpha=0.05$）?

分析：

电子元件的寿命总体均值 $μ$ 未知，故可以采用 $T$ 检验。

需要检验元件的平均寿命是否大于225小时（$\mu>225$），即右侧检验。

$H_{0}:\mu\leq225$ $H_{1}:\mu>225$

$$T=\frac{\overline{x}-\mu}{s/\sqrt{n}}$$

需要自己计算出寿命 $x$ 的样本均值、样本标准差。

实现代码:

import math
import numpy as np
from scipy import stats


def t_test(data, n, test_time):
    # 计算样本均值和标准差
    sample_mean = np.mean(data)
    sample_std = np.std(data, ddof=1)

    # 计算标准误差和 T 统计量
    SE = sample_std / math.sqrt(n)
    t_score = (sample_mean - test_time) / SE

    # 计算临界值
    alpha = 0.05
    df = n - 1
    critical_value = stats.t.ppf(1 - alpha, df)

    # 计算 p 值
    p_value = stats.t.sf(np.abs(t_score), df)

    return t_score, critical_value, p_value


if __name__ == "__main__":
    # 样本数据
    data = [159, 280, 101, 212, 224, 379, 179, 264, 222, 362, 168, 250, 149, 260, 485, 170]
    # 总体和样本参数
    test_time = 225  # 检验时间
    n = 16  # 样本量
    alpha = 0.05

    t_score, critical_value, p_value = t_test(data, n, test_time)

    print(f"T值: {t_score:.4f}")
    print(f"临界值: {critical_value:.4f}")
    print(f"P值: {p_value:.4f}")

    print()
    print("从临界值角度:")
    if t_score < critical_value:
        print("落在非拒绝域,接受原假设，元件的平均寿命 ≤ 225小时")
    else:
        print("落在拒绝域,拒绝原假设，元件的平均寿命 > 225小时")

    print("从P值角度:")
    if p_value < alpha:
        print("拒绝原假设，元件的平均寿命 > 225小时")
    else:
        print("接受原假设，元件的平均寿命 ≤ 225小时")

输出结果：

T值: 0.6685
临界值: 1.7531
P值: 0.2570

从临界值角度:
落在非拒绝域,接受原假设，元件的平均寿命 ≤ 225小时
从P值角度:
接受原假设，元件的平均寿命 ≤ 225小时

经过$T$ 检验后，

从临界值角度：

$$t=\frac{\overline{x}-\mu_{0}}{s/\sqrt{n}}=0.669<t_{0.05}(16)=1.753$$

t值小于自由度为16的 $T$ 分布的上0.05分位点，落在非拒绝域，故接受原假设。

从P值角度：

$$P_{\mathrm{value}}=P(T>0.669)=0.257>\alpha=0.05$$

如果拒接原假设，犯弃真错误的概率达0.257，大于 $\alpha=0.05$ ，故接受原假设。

综上，接受原假设，元件的平均寿命 ≤ 225小时。

实验五

描述：给出两个.csv文件，数据分别为某两股票为期一年的日交易信息

请根据每只股票的收盘价与开盘价计算相应每天的波动状况，计算一年来那只股票的每天差额波动大？

每天差额值的分布是否服从正态分布，是否是对称的？

分析：

数据的离散程度：

要比较股票的差额波动程度，可以计算差额的方差和变异系数。

方差：

$$\sigma^2=\frac{\sum(X-\mu)^2}N$$

标准差：

$$\sigma=\sqrt{\sigma}$$

平均值：

$$\mu=\frac{\sum_{i=1}^{n}x_{i}}{n}$$

变异系数：

$$c_v=\frac\sigma\mu$$

正态分布检验：

通过 $cipy.stats.normaltest$ 计算 $p$ 值

通过 $.skew()$ 计算 $symmetry$

$$p\begin{cases}<0.05\text{ 不服从正态分布}\\\\\geq0.05\text{ 服从正态分布}\end{cases}$$

$$\text{symmetry}\begin{cases}>0&\text{右偏}\\=0&\text{对称}\\<0&\text{左偏}\end{cases}$$

实现代码:

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from scipy.stats import normaltest, norm


def stock_analysis(path, num):
    # 读取 CSV 文件
    df = pd.read_csv(path, encoding='GBK')
    # 计算收盘价与开盘价的差额
    df['差额'] = df['收盘价'] - df['开盘价']

    # 方差
    var = df['差额'].var()
    # 标准差
    std = df['差额'].std()
    # 平均值
    mean = df['差额'].mean()
    # 变异系数
    cv = (std - mean) * 100 if mean else 0

    # 检验差额是否服从正态分布
    stat, p_value = normaltest(df['差额'])
    # 检验差额是否对称
    symmetry = df['差额'].skew()

    print(f"股票{num}的波动情况：")
    print(f"方差: {var:.4f}")
    print(f"变异系数: {cv:.4f}")
    print()

    # 绘制正态分布图
    plt.figure()
    plt.hist(df['差额'], bins=30, density=True, alpha=0.6, color='g')
    xmin, xmax = plt.xlim()
    x = np.linspace(xmin, xmax, 100)
    p = norm.pdf(x, mean, std)
    plt.plot(x, p, 'k', linewidth=2)
    title = f'股票{num}的差额正态分布图'
    plt.title(title)
    plt.show()

    return var, cv, p_value, symmetry


def normal_analysis(p_value, symmetry, num):
    if p_value < 0.05:
        print(f"股票{num}：差额不服从正态分布")
    else:
        print(f"股票{num}：差额服从正态分布")

    if symmetry == 0:
        print(f"股票{num}：差额对称")
    elif symmetry > 0:
        print(f"股票{num}：差额右偏")
    else:
        print(f"股票{num}：差额左偏")
    print()


if __name__ == "__main__":

    # 加载中文字体
    plt.rcParams['font.sans-serif'] = ['SimSun']
    # 解决正负号问题
    plt.rcParams["axes.unicode_minus"] = False

    path1 = 'daily_20240428200156.csv'
    var1, cv1, p_value1, symmetry1 = stock_analysis(path1, 1)
    path2 = 'daily_20240428200125.csv'
    var2, cv2, p_value2, symmetry2 = stock_analysis(path2, 2)

    print("结论：", end="")
    if var1 > var2:
        print("股票1每天差额波动较大\n")
    elif var1 < var2:
        print("股票2每天差额波动较大\n")
    else:
        print("两只股票每天差额波动相同\n")

    print("正态分析：")
    normal_analysis(p_value1, symmetry1, 1)
    normal_analysis(p_value2, symmetry2, 2)

输出结果：

股票1的波动情况：
方差: 0.0670
变异系数: 22.6865

股票2的波动情况：
方差: 5.0120
变异系数: 202.7938

结论：股票2每天差额波动较大

正态分析：
股票1：差额不服从正态分布
股票1：差额右偏

股票2：差额不服从正态分布
股票2：差额右偏

最后画出两只股票的正态分布图：

实验六

描述：

根据方差分析的描述，编写一个单因素方差分析的函数。

输入为某二维数据（列代表因素各水平值）

输出为检验统计量 $F$ 的值

分析：

0.数据结构

$A_1$	$A_2$	···	$A_k$
$x_{11}$	$x_{21}$	···	$x_{k1}$
$x_{12}$	$x_{12}$	···	$x_{k2}$
⋮	⋮	⋮	⋮
$x_{1n_1}$	$x_{2n_2}$	···	$x_{kn_k}$

​ 在单因素方差分析中，用 $A$ 表示因素，因素的 $k$ 个水平 (总体)分别用 $A_1, A_2, ···, A_k$ 表示，每个观测值用 $x_{ij}$ ( $i=1, 2, ···, k; j=1, 2, ···, n_i$ )表示，即 $x_{ij}$ 表示第 $i$ 个水平(总体)的第 $j$ 个观测值。值得注意的是，每个因素的观测值数量可能不同，即 $n_i$ 不同。

1.提出假设

$H_{0}: \mu_{1}=\mu_{2}=\cdots=\mu_{i}=\cdots=\mu_{k}$ 自变量对因变量没有显著影响

$H_{1}: \mu_{i} (i=1, 2, \cdots, k)$ 自变量对因变量有显著影响

式中，$\mu_{i}$为第 $i$ 个总体的均值。

2.构造检验的统计量

(1) 计算各样本的均值

​ 假定从第 $i$ 个总体中抽取一个容量为 $n_i$ 的简单随机样本，令 $\overline{x}_{i}$ 为第 $i$ 个总体的样本均值，则有

$$\\\bar{x}_i=\frac{\sum_{j=1}^{n_i}x_{ij}}{n_i},\quad i=1,2,\cdots,k$$

式中， $n_i$ 为第 $i$ 个总体的样本量； $x_{ij}$ 为第 $i$ 个总体的第 $j$ 个观测值。

(2) 计算全部观测值的总均值

​ 全部观测值的总和除以观测值总个数的结果。令总均值为 $\overline{\overline{x}}$ ，则有

$$\overline{\overline{x}}=\frac{\sum_{i=1}^k\sum_{j=1}^{n_i}x_{ij}}n=\frac{\sum_{i=1}^kn_i\overline{x}_i}n$$

式中， $n=n_1+n_2+\cdotp\cdotp\cdotp+n_k$

(3) 计算各误差平方和

​ 分别计算三个误差平方和，它们是总平方和、组间平方和（因素平方和）、组内平方和（误差平方和或残差平方和）。

1. 总平方和

记为 $SST$ ，它是全部观测值 $x_{ij}$ 与总均值 $\overline{x}$ 的误差平方和，反映各所有观测值与总均值之间的差异

$$\text{SS}T=\sum_{i=1}^{k}\sum_{j=1}^{n_{i}}{(x_{ij}-\overline{x})^{2}}$$

2. 组间平方和

记为 $SSA$ ，它是各组均值 $\overline{x}_{i}$ ( $i=1, 2, ···, k$ ) 与总均值 $\overline{x}$ 的误差平方和，反映各样本均值之间的差异程度

$$SSA = \sum_{i=1}^{k}n_{i}(\overline{x}_{i}-\overline{x})^{2}$$

3. 组内平方和

记为 $SSE$ ，它是每个水平或组的各样本数据与其组均值的误差平方和，反映每个样本各观测值的离散状况

$$\mathrm{SS}E=\sum_{i=1}^{k}\sum_{i=1}^{n_{i}}(x_{ij}-\overline{x}_{i})^{2}$$

4. 计算统计量

为了消除观测值多少对误差平方和大小的影响，需要将其平均，也就是用各平方和除以它们所对应的自由度，这一结果称为均方。

$SST$ 的自由度为 $n-1$ ，其中 $n$ 为全部观测值的个数。

$SSA$ 的自由度为 $k-1$ ，其中 $k$ 为因素水平（总体）的个数。

$SSE$ 的自由度为 $n-k$ 。

由于要比较的是组间均方和组内均方之间的差异，所以通常只计算 $SSA$ 的均方和 $SSE$ 的均方。

$SSA$ 的均方也称组间均方或组间方差，记为 $MSA$

$$M\text{SA}=\frac{\text{组间平方和}}{\text{自由度}} = \frac { \mathrm{SS}A}{k-1}$$

$SSE$ 的均方也称为组内均方或组内方差，记为 $MSE$

$$MSE=\frac{\text{组内平方和}}{\text{自由度}} = \frac { S S E }{ n - k}$$

将上述 $MSA$ 和 $MSE$ 进行对比，即得到所需要的检验统计量 $F$ 。

当 $H_0$ 为真时，二者的比值服从分子自由度为 $k-1$ 、分母自由度为 $n-k$ 的 $F$ 分布，即

$$F=\frac{MSA}{MSE}\sim F(k-1,n-k)$$

3.作出统计决策

​ 正常的单因素方差分析还需根据 $F$ 分布表，作出对原假设 $H_0$ 的决策，此处只需求出统计量 $F$ 的值即可，不过多赘述。

实现代码:

import numpy as np


def one_way_analysis_of_variance(data):
    # 初始化组间平方和(SSA)和组内平方和(SSE)
    SSA = 0
    SSE = 0

    # 计算每组非零数据的数量
    nonzero_group_num = np.count_nonzero(data != 0, axis=0)
    # 计算非零数据的总数
    nonzero_total_num = data.size - np.count_nonzero(data == 0)
    # 计算非零数据的总均值
    nonzero_total_mean = np.mean(data[data != 0])

    # 对每一列（每个组）进行操作
    for i, group_data in enumerate(data.T):
        # 排除0值后计算当前组的均值
        group_nonzero_data = group_data[group_data != 0]
        group_mean = np.mean(group_nonzero_data)

        # 更新组间平方和(SSA)和组内平方和(SSE)
        if group_nonzero_data.size > 0:  # 确保组内至少有一个非零观测值
            SSA += group_nonzero_data.size * (group_mean - nonzero_total_mean) ** 2
            SSE += np.sum((group_nonzero_data - group_mean) ** 2)

    # 计算自由度
    SSA_freedon_degree = data.shape[1] - 1  # SSA自由度，列数-1
    SSE_freedon_degree = nonzero_total_num - data.shape[1]  # SSE自由度，非0数据总数-列数

    # 计算均方(MSA和MSE)
    MSA = SSA / SSA_freedon_degree
    MSE = SSE / SSE_freedon_degree

    # 计算F统计量
    F_value = MSA / MSE

    return F_value


# 示例数据
data = np.array([  # 0代表无数据
    [57, 68, 31, 44],
    [66, 39, 49, 51],
    [49, 29, 21, 65],
    [40, 45, 34, 77],
    [34, 56, 40, 58],
    [53, 51, 0, 0],
    [44, 0, 0, 0]
])

if __name__ == "__main__":
    F_value = one_way_analysis_of_variance(data)
    print(f"检验统计量 F 的值为: {F_value:.6f}")

输出结果：

检验统计量 F 的值为: 3.406643

实验七

描述：

编写两个 一元线性回归中的函数

第一个函数求相关系数

第二个函数求一次项系数

分析：

输入两个一维向量序列，分别求出相关系数和一次项系数

1.相关系数

相关系数是研究变量之间线性相关程度的量。

$$r\left(X,Y\right)=\frac{Cov\left(X,Y\right)}{\sqrt{Var\left[X\right]Var\left[Y\right]}}$$

其中 $Cov\left(X,Y\right)$ 是协方差， $Var$ 是方差。

2.一次项系数

参数的最小二乘估计

$$\hat{y}_{i}=\hat{\beta}_{0}+\hat{\beta}_{1}x_{i}$$

$$\hat{\beta}_{1}=\frac{n\sum_{i=1}^{n}x_{i}y_{i}-\sum_{i=1}^{n}x_{i}\sum_{i=1}^{n}y_{i}}{n\sum_{i=1}^{n}x_{i}^{2}-\left(\sum_{i=1}^{n}x_{i}\right)^{2}}$$

其中 $\hat{\beta}_{1}$ 即为一次项系数

实现代码:

def correlation_coefficient(x, y):
    # 计算均值
    mean_x = sum(x) / len(x)
    mean_y = sum(y) / len(y)

    # 计算协方差
    covariance = sum((x[i] - mean_x) * (y[i] - mean_y) for i in range(len(x)))

    # 计算分母
    var_x = sum((x[i] - mean_x) ** 2 for i in range(len(x)))
    var_y = sum((y[i] - mean_y) ** 2 for i in range(len(y)))
    denominator = (var_x * var_y) ** 0.5

    # 计算相关系数
    r = covariance / denominator
    return r


def linear_regression(x, y):
    sum_x = sum(x)
    sum_y = sum(y)
    sum_xy = sum(x[i] * y[i] for i in range(len(x)))
    sum_x_square = sum(x[i] ** 2 for i in range(len(x)))

    # 计算一次项系数
    a = (len(x) * sum_xy - sum_x * sum_y) / (len(x) * sum_x_square - sum_x ** 2)

    return a


# 测试数据
x = [67.3, 111.3, 173.0, 80.8, 199.7, 16.2, 107.4, 185.4, 96.1, 72.8, 64.2, 132.2, 58.6, 174.6, 263.5, 79.3, 14.8, 73.5,
     24.7, 139.4, 368.2, 95.7, 109.6, 196.2, 102.2]
y = [0.9, 1.1, 4.8, 3.2, 7.8, 2.7, 1.6, 12.5, 1.0, 2.6, 0.3, 4.0, 0.8, 3.5, 10.2, 3.0, 0.2, 0.4, 1.0, 6.8, 11.6, 1.6,
     1.2, 7.2, 3.2]

if __name__ == "__main__":
    print(f"相关系数为：{correlation_coefficient(x, y):.6f}")
    print(f"一次项系数为：{linear_regression(x, y):.6f}")

输出结果：

相关系数为：0.843571
一次项系数为：0.037895