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'
