感谢Dr.fish的耐心讲解和细致回答。
本次课的课后作业如下：
分别用 T 分布 和 bootsrrap 方法求年均降水量数据在置信度为95%的置信区间
上代码
# 导入分析包及数据
import scipy.stats
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%config InlineBackend.figure_format = 'retina'
df = pd.read_csv('rainfall.csv', header=None)
df_size = df.iloc[:,0] # 取全部数据
做作业
Q1. 利用 T 分布求年均降水量数据在置信度为95%的置信区间
通过两种方式进行T分布计算
# 直接计算T分布，置信度95%下年均降水量的置信区间
df_std = df.std() # 计算样本标准差
sample_mean = df_size.mean() # 计算样本均值
sample_size = len(df_size)
t_score = scipy.stats.t.pdf(0.025 , sample_size - 1)
margin_error = t_score * df_std / np.sqrt(sample_size)
lower_limit = sample_mean - margin_error
upper_limit = sample_mean + margin_error
print '95%% Confidence Interval: ( %.1f, %.1f)' % (lower_limit, upper_limit)
#计算结果
95% Confidence Interval: ( 795.9, 804.3)
# 定义函数计算T分布，置信度95%下年均降水量的置信区间
def ci_t(data, df_std, confidence):
sample_mean = np.mean(data)
sample_size = len(data)
alpha = (1 - confidence) / 2
t_score = scipy.stats.t.pdf(alpha , sample_size - 1)
ME = t_score * df_std / np.sqrt(sample_size)
lower_limit = sample_mean - ME
upper_limit = sample_mean + ME
return (lower_limit , upper_limit)
#输入变量
ci_t(df_size, df_std, 0.95)
#计算结果
(795.89206, 804.291398)
Q2. 利用 bootsrrap 分布求年均降水量数据在置信度为95%的置信区间
#定义计算均值和bootstrap函数
def bootstrap_mean(data):
return np.mean(np.random.choice(data, size = 10)) # 从数据data中有放回抽样，每次抽10个，并返回样本均值
def draw_bootstrap(data, times = 1):
bs_mean = np.empty(times) #初始化长度为times的空数组
for i in range(times): #进行times次抽样，将每次得到的样本均值存储在bs_mean中
bs_mean[i] = bootstrap_mean(data)
return bs_mean
#输入变量
bs_mean = draw_bootstrap(df_size, 10000) #执行有放回抽样一万次
#计算置信区间
np.percentile(bs_mean, [2.5, 97.5])
#计算结果
array([ 725.878, 876.481])
出个图看一下
#计算 bin size
IQR = df_size.quantile(0.75) - df_size.quantile(0.25)
bin_size = int(2 * IQR / len(df_size)**(1/ 3))
#绘图
plt.hist(bs_mean, bins = bin_size, normed = True, color = 'r', alpha = 0.5, rwidth = 0.9, label='T distribution')
plt.axvline(np.percentile(bs_mean, 2.5),color = 'g', alpha = 0.5)
plt.axvline(np.percentile(bs_mean, 97.5),color = 'g', alpha = 0.5)
plt.legend() #显示数据标签
plt.show()
bootstrap
关于bin_size
之前一直纠结 histogram 到底取多少 bins 合适，少了会把异常掩埋，多了……（暂时没想到问题在哪，尬笑）。后来有次问 Dr.fish，她分享了 wiki 给我。作为一个丢三落四的技术白，如果不存起来肯定就找不到了，所以丢到作业里，照着算一下再存一下，就算下次用忘了也知道去哪儿考古，吼吼（小白猜你们一定能体会我此刻的心情~）
bin_size
