商业分析第五次课作业-0808（求置信度95%的降水概率）

感谢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