import warnings
warnings.filterwarnings('ignore')
import pandas as pd 
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

均数与标准差

某省的高考分数经过标准化以后，最低分为100分，最高分为900分，平均分为500分，标准差为100分。用计算机模拟从该总体中随机抽取20名考生的分数见下表。试计算考生分数的均数与标准差。

data = {
'考生号':[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20],
 '分数': [456, 594, 611, 336, 298, 394, 464, 336, 513, 553, 541, 478, 306, 516, 456, 452, 431, 531, 435, 552]
}
df = pd.DataFrame(data)

mean = df['分数'].mean()
std_dev = df['分数'].std()

print("均值:", mean)
print("标准差:", std_dev)
均值: 462.65
标准差: 92.4082902409906

统计描述

从幼儿园大班随机抽取12名6周岁女童，测得身高（cm）。试进行统计描述。

data = {
    '编号': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12],
    '身高（cm）': [125.2, 135.3, 122.9, 131.6, 121.1, 141.5, 132.1, 112.8, 104.6, 131.2, 125.9, 126.1]
}

df = pd.DataFrame(data)
statistics = df['身高（cm）'].describe()
print(statistics)
count     12.000000
mean     125.858333
std        9.947997
min      104.600000
25%      122.450000
50%      126.000000
75%      131.725000
max      141.500000
Name: 身高（cm）, dtype: float64

统计描述提供了以下信息：

count：非缺失值的个数。 mean：平均值。 std：标准差。 min：最小值。 25%：第一四分位数，即将数据从小到大排序后，处于中间位置的值。 50%：中位数，即将数据从小到大排序后，处于中间位置的值。 75%：第三四分位数，即将数据从小到大排序后，处于中间位置的值。 max：最大值。

总体均数估计

某县1998年抽样调查了500户农民家庭的年化纤布消费量，得到均数为3.55米，标准差为1.03米。试估计该县1998年农民家庭年化纤布消费量的总体均数。 由于样本数据的均数为3.55米，标准差为1.03米，我们可以认为样本数据的分布是正态分布。因此，估计总体均数的值就是样本均数3.55米

# 样本均数
mean = 3.55
# 样本标准差
std_dev = 1.03
# 总体均数的估计值等于样本均数
estimate = mean
print("总体均数的估计值为：", estimate)
总体均数的估计值为： 3.55

总体率估计

为了解某地新生儿畸形的发生率，某单位调查了该地3009名活产新生儿，诊断出畸形者29名，占0.96%。试估计该地活产新生儿的畸形率。该地活产新生儿中诊断出畸形者29名，占0.96%。我们可以使用二项分布来估计该地活产新生儿的畸形率。二项分布的概率公式为：

P(X=k) = C(n, k) * p^k * (1-p)^(n-k)

在这个问题中，n=3009，k=29，p=0.96%。我们需要计算P(X>=k)，即活产新生儿中至少有k个畸形者的概率。这个概率的值等于1减去活产新生儿中畸形数小于k的概率。即：

P(X>=k) = 1 - P(X<k) = 1 - ∑[C(n, i) * p^i * (1-p)^(n-i)]，其中i=0,1,...,k-1。

import scipy.stats as stats
# 活产新生儿总数
n = 3009
# 诊断出畸形的新生儿数
k = 29
# 畸形率
p = 0.96 / 100
# 计算P(X>=k)的值
prob = 1 - stats.binom.cdf(k - 1, n, p)
print("该地活产新生儿的畸形率估计值为：", prob * 100)
该地活产新生儿的畸形率估计值为： 51.660804398075754

样本均数与总体均数的比较

据大量调查知，健康成年男子脉搏的均数为72次/分，某医生在山区随机调查了25名健康成年男子，其脉搏均数为75.5次/分，标准差为6.5次/分，能否认为该山区成年男子的脉搏高于一般人群？

import numpy as np
from scipy import stats

# 假设检验
population_mean = 72  # 总体均数
sample_mean = 75.5  # 样本均数
sample_std = 6.5  # 样本标准差
n = 25  # 样本容量

# 计算t统计量
t_statistic = (sample_mean - population_mean) / (sample_std / np.sqrt(n))

# 计算p值
p_value = 1 - stats.t.cdf(t_statistic, df=n-1)

# 输出结果
if p_value < 0.05:
    print("根据t检验，可以认为该山区成年男子的脉搏高于一般人群，p值为", p_value)
else:
    print("根据t检验，无法认为该山区成年男子的脉搏高于一般人群，p值为", p_value)
根据t检验，可以认为该山区成年男子的脉搏高于一般人群，p值为 0.006364545798554722

配对设计的两样本均数的比较

欲研究某药物对血红蛋白含量是否有影响，观察了9例患者治疗前后血红蛋白的变化，数据如下表。试问，该药物治疗前后血红蛋白含量有无变化？

import numpy as np
from scipy import stats

# 原始数据
data = np.array([[122, 145], [113, 128], [141, 156], [123, 122], [105, 121],
                 [124, 105], [144, 123], [115, 101], [117, 127]])

# 计算差值
diff = data[:, 1] - data[:, 0]

# 进行假设检验
t_statistic, p_value = stats.ttest_1samp(diff, 0)

# 输出结果
if p_value < 0.05:
    print("根据t检验，药物治疗前后血红蛋白含量有显著变化，p值为", p_value)
else:
    print("根据t检验，药物治疗前后血红蛋白含量没有显著变化，p值为", p_value)
根据t检验，药物治疗前后血红蛋白含量没有显著变化，p值为 0.6475376620106789

两个样本均数比较（成组设计）

欲研究某药物对血红蛋白含量是否有影响，把18例患者随机分为实验组（用该药物治疗）和对照组（用对血红蛋白无影响的标准药物治疗），每组各9例，治疗后两组患者血红蛋白含量如下表。试问，该药物是否影响血红蛋白含量？ 实验组 122 113 141 123 105 124 144 115 117 对照组 148 129 156 122 121 105 123 100 126

import numpy as np
from scipy import stats

# 实验组数据
experimental_group = [122, 113, 141, 123, 105, 124, 144, 115, 117]
# 对照组数据
control_group = [148, 129, 156, 122, 121, 105, 123, 100, 126]

# 假设检验
t_statistic, p_value = stats.ttest_ind(experimental_group, control_group)

# 输出结果
if p_value < 0.05:
    print("根据独立两样本t检验，可以认为该药物对血红蛋白含量有影响，p值为", p_value)
else:
    print("根据独立两样本t检验，无法认为该药物对血红蛋白含量有影响，p值为", p_value)
根据独立两样本t检验，无法认为该药物对血红蛋白含量有影响，p值为 0.697947145973673

多个样本均数比较（成组设计）

欲研究药物A、B对血红蛋白含量是否有影响，把15例患者随机分为A药组（用A药物治疗）、B药组（用B药物治疗）和对照组（用安慰剂治疗），每组各5例，治疗后各组患者血红蛋白含量如下表。试问，药物A、B是否影响血红蛋白含量？ A药组 122 113 141 123 105 B药组 144 126 156 122 121 对照组 101 111 113 100 101

import numpy as np
from scipy import stats

# A药组数据
group_a = [122, 113, 141, 123, 105]
# B药组数据
group_b = [144, 126, 156, 122, 121]
# 对照组数据
control_group = [101, 111, 113, 100, 101]

# 单因素方差分析
f_statistic, p_value = stats.f_oneway(group_a, group_b, control_group)

# 输出结果
if p_value < 0.05:
    print("根据单因素方差分析，可以认为药物A、B对血红蛋白含量有影响，p值为", p_value)
else:
    print("根据单因素方差分析，无法认为药物A、B对血红蛋白含量有影响，p值为", p_value)
根据单因素方差分析，可以认为药物A、B对血红蛋白含量有影响，p值为 0.011236776709940352

配伍组设计多个样本均数比较

为研究药物A、B对血红蛋白含量是否有影响，把15例患者根据性别、年龄、文化程度等因素分为5个区组，即每个区组的3个人性别相同、年龄和文化程度相近，再把每个区组的3个人随机分配到A药组（用A药物治疗）、B药组（用B药物治疗）和对照组（用安慰剂治疗）中。治疗后各组患者血红蛋白含量如下表。试问，药物A、B是否影响血红蛋白含量？ 区组号 1 2 3 4 5 A药组 1 122 113 141 123 105 B药组 2 144 126 156 122 121 对照组3 103 110 115 100 101

import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols

# 创建数据框
data = pd.DataFrame({
    'Group': ['A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'C', 'C', 'C', 'C', 'C'],
    'Value': [122, 113, 141, 123, 105, 144, 126, 156, 122, 121, 103, 110, 115, 100, 101]
})

# 执行单因素方差分析
model = ols('Value ~ Group', data=data).fit()
anova_table = sm.stats.anova_lm(model)

# 打印方差分析结果
print(anova_table)
df       sum_sq     mean_sq         F   PR(>F)
Group      2.0  1963.333333  981.666667  6.359318  0.01309
Residual  12.0  1852.400000  154.366667       NaN      NaN

样本率与总体率的比较

据大量调查知，一般溃疡病患者中有20%发生胃出血症状，某医生观察245例70岁以上溃疡病人，其中75例发生出血症状，问老年患者与一般患者胃出血发生率是否不同？

import statsmodels.api as sm

# 定义总样本数和发生出血症状的样本数
n1 = 245  # 总样本数
x1 = 75   # 发生出血症状的样本数

# 定义比例
p0 = 0.2  # 一般溃疡病患者中的出血发生率

# 执行比例检验
result = sm.stats.proportions_ztest(x1, n1, p0)

# 打印检验结果
print("Z-score:", result[0])
print("P-value:", result[1])
Z-score: 3.6041370564350945
P-value: 0.0003131917647302031

完全随机设计两个样本率的比较（四格表资料）

为研究甲乙两种药物对胃溃疡的治疗效果，选择了128名病例，随机分为两组，治疗结果结果如表1。问甲乙两种药物对胃溃疡的疗效有无差别？ 表1 甲乙两种药物对胃溃疡的疗效 组别 结果 结果 合计 痊愈 无效
A药物 60 4 64 B药物 48 16 64 合计 108 20 128

import numpy as np
import scipy.stats as stats

# 创建观察矩阵
observed = np.array([[60, 4],
                     [48, 16]])

# 执行卡方检验
chi2_stat, p_val, _, _ = stats.chi2_contingency(observed)

# 打印检验结果
print("Chi-square statistic:", chi2_stat)
print("P-value:", p_val)
Chi-square statistic: 7.170370370370371
P-value: 0.007411747131009725

多个样本率的比较

为研究某药物治疗胃溃疡的疗效，把105名患者随机分为三组，得资料如表1，问不同剂量的疗效是否相同？ 表1 三种不同剂量的治疗结果 剂量 有效 无效 合计 小剂量 19 11 30 中剂量 41 9 50 大剂量 24 1 25 合计 82 23 105

import numpy as np
import scipy.stats as stats

# 创建观察矩阵
observed = np.array([[19, 11],
                     [41, 9],
                     [24, 1]])

# 执行卡方检验
chi2_stat, p_val, _, _ = stats.chi2_contingency(observed)

# 打印检验结果
print("Chi-square statistic:", chi2_stat)
print("P-value:", p_val)
Chi-square statistic: 9.333333333333334
P-value: 0.009403562551495204

单向有序分类资料的假设检验

某研究得资料如表1，问2种药物的疗效是否相同？ 表1 2种药物疗效的观察结果 药物 疗效 合计 治愈1 显效2 好转3 无效4 A药物1 26 23 10 1 60 B药物2 12 15 21 12 60 合计 33 37 31 13 120

import numpy as np
import scipy.stats as stats

# 创建观察矩阵
observed = np.array([[26, 23, 10, 1],
                    [12, 15, 21, 12]])

# 执行卡方检验
chi2_stat, p_val, _, _ = stats.chi2_contingency(observed)

# 打印检验结果
print("Chi-square statistic:", chi2_stat)
print("P-value:", p_val)
Chi-square statistic: 20.053023377301816
P-value: 0.00016550122567394294

某研究得资料如表2，问病型与患者痰液中SB的含量是否有关系？ 表2 病型与患者痰液中SB含量的关系 病型 SB含量 合计 - + ++ +++ A型 12 22 34 22 90 B型 11 12 32 2 57 C型 10 34 23 1 68 D型 5 2 23 3 33

import numpy as np
import scipy.stats as stats

# 创建观察矩阵
observed = np.array([[12, 22, 34, 22],
                     [11, 12, 32, 2],
                     [10, 34, 23, 1],
                     [5, 2, 23, 3]])

# 执行卡方检验
chi2_stat, p_val, _, _ = stats.chi2_contingency(observed)

# 打印检验结果
print("Chi-square statistic:", chi2_stat)
print("P-value:", p_val)
Chi-square statistic: 51.24036530603176
P-value: 6.290216614589131e-08

相关分析

为了解城市儿童年龄与身高的关系，在某小学随机抽取8名6—12岁儿童，测得身高如下表。问儿童身高与年龄之间是否相关？ 编号 1 2 3 4 5 6 7 8 年龄（岁）x 6.2 7.0 10.2 11.0 12.1 9.5 8.2 6.5 身高（cm）y 135 139 143 150 155 141 140 137

import numpy as np
import scipy.stats as stats

# 创建年龄和身高的数组
age = np.array([6.2, 7.0, 10.2, 11.0, 12.1, 9.5, 8.2, 6.5])
height = np.array([135, 139, 143, 150, 155, 141, 140, 137])

# 计算Pearson相关系数
r, p = stats.pearsonr(age, height)

# 打印结果
print("Pearson correlation coefficient:", r)
print("P-value:", p)
Pearson correlation coefficient: 0.9357981688909468
P-value: 0.0006301328906580182

回归分析

为了解城市儿童年龄与身高的关系，在某小学随机抽取8名6—12岁儿童，测得身高如下表。试建立直线回归方程。 编号 1 2 3 4 5 6 7 8 年龄（岁） 6.2 7.0 10.2 11.0 12.1 9.5 8.2 6.5 身高（cm） 135 139 143 150 155 141 140 137

import numpy as np
import statsmodels.api as sm

# 创建年龄和身高的数组
age = np.array([6.2, 7.0, 10.2, 11.0, 12.1, 9.5, 8.2, 6.5])
height = np.array([135, 139, 143, 150, 155, 141, 140, 137])

# 添加常数项（截距）
age_with_const = sm.add_constant(age)

# 拟合线性回归模型
model = sm.OLS(height, age_with_const).fit()

# 打印回归系数
print(model.summary())

# 输出回归方程
print("regression equation: y = {:.2f} + {:.2f}x".format(model.params[0], model.params[1]))
OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.876
Model:                            OLS   Adj. R-squared:                  0.855
Method:                 Least Squares   F-statistic:                     42.28
Date:                Sat, 23 Dec 2023   Prob (F-statistic):           0.000630
Time:                        22:31:36   Log-Likelihood:                -17.766
No. Observations:                   8   AIC:                             39.53
Df Residuals:                       6   BIC:                             39.69
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const        117.0605      4.017     29.141      0.000     107.231     126.890
x1             2.8786      0.443      6.502      0.001       1.795       3.962
==============================================================================
Omnibus:                        0.835   Durbin-Watson:                   2.741
Prob(Omnibus):                  0.659   Jarque-Bera (JB):                0.659
Skew:                          -0.454   Prob(JB):                        0.719
Kurtosis:                       1.925   Cond. No.                         40.5
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
regression equation: y = 117.06 + 2.88x


d:\programdata\anaconda3\lib\site-packages\scipy\stats\stats.py:1394: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=8
  "anyway, n=%i" % int(n))

生存分析

import numpy as np
from lifelines import KaplanMeierFitter

# 创建数据数组
a = np.array([200, 155, np.inf, np.inf, 300, 166, 215, 312, np.inf, 213, 166, 222])
b = np.array([135, 139, np.inf, 143, 150, 155, 141, np.inf, 140, 137])

# 创建事件数组和状态数组
events_a = np.where(a == np.inf, 0, 1)
events_b = np.where(b == np.inf, 0, 1)
statuses_a = np.where(a == np.inf, 0, a)
statuses_b = np.where(b == np.inf, 0, b)

# 创建Kaplan-Meier对象并拟合数据
kmf = KaplanMeierFitter()
kmf.fit(statuses_a, events=events_a, label='A方案')
kmf.fit(statuses_b, events=events_b, label='B方案')

# 绘制Kaplan-Meier曲线
kmf.plot()
---------------------------------------------------------------------------

ModuleNotFoundError                       Traceback (most recent call last)

~\AppData\Local\Temp\ipykernel_17728\2082065611.py in <module>
      1 import numpy as np
----> 2 from lifelines import KaplanMeierFitter
      3 
      4 # 创建数据数组
      5 a = np.array([200, 155, np.inf, np.inf, 300, 166, 215, 312, np.inf, 213, 166, 222])


ModuleNotFoundError: No module named 'lifelines'