实践-python简单实现参数估计
作业：https://mp.weixin.qq.com/s/8egc4QE6MmME0AS4FKSlOg
下面动手实践 主要是借鉴：
https://blog.csdn.net/polarislove36/article/details/78922045?utm_source=blogxgwz4
https://blog.csdn.net/maplepiece1999/article/details/103676556
import pandasas pd
import numpyas np
import matplotlib.pyplotas plt
from scipyimport stats
import scipy
#区间估计
#区间估计能够提供有关估计精度的信息。
#我们通过点估计值加上和减去一个边际误差（margin of error）的值，可以构建出一个总体参数的区间估计
#8.1大样本估计
#假定σσ 已知的大样本情况（n>=30）：
#总计均值的区间估计
#x¯¯¯±zα/2σn−−√
#x¯±zα/2σn
#其中1-αα 是置信系数，zα/2zα/2 是当标准正态分布的上侧面积为 α/2α/2 时的z的值
#这种情况需要知道总体的标准差，如果不知道的情况需要用总体标准差的样本点估计值s
#如果是小样本，并且知道总体的标准差σσ 需要假设总体满足正态分布，方法和大样本相同
# #如果是小样本，不知道总体的标准差σσ，只知道样本的标准差s需要假设总体满足正态分布，并且使用t分布
#t分布是一系列类似的概率分布，依赖于自由度，自由度越大，越接近正态分布
#小样本且假设总体满足正态分布：用s估计σσ 的小样本情况（n<30）:
#x¯¯¯±tα/2sn−−√
# 总体 方差的区间估计
#df2 = df.groupby(['Embarked'])
#df2.describe()
# = pd.read_csv('rainfall.csv',header = None)
data = pd.read_excel('D:/Users/Desktop/Desktop/作业/data.xlsx', usecols=[1])
data = data.iloc[:,0]
print(data.mean())
# 抽取100个样本
data_sam = data.sample(100)
x = data_sam.mean()#求抽样平均值
data_sam.describe()
#计算置信区间
#计算标准差
#pandas.std() 默认是除以n-1,无偏；加上参数ddof=0，即为有偏；DataFrame的describe()中就包含有std()；
#numpy.std() 求标准差的时候默认是除以 n 的，即是有偏的，np.std无偏样本标准差方式为加入参数ddof = 1
# 正态分布下的置信区间
def norm_conf (data,confidence=0.95):
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html
sample_mean = np.mean(data)
sample_std = np.std(data,ddof=1)
sample_size =len(data)
conf_intveral = scipy.stats.norm.interval(confidence, loc=sample_mean, scale=sample_std)
print(conf_intveral)
return conf_intveral
# T分布下的置信区间
def ttest_conf(data, confidence=0.95):
sample_mean = np.mean(data)
sample_std = np.std(data, ddof=1)
sample_size =len(data)
conf_intveral = scipy.stats.t.interval(confidence, df=(sample_size -1), loc=sample_mean, scale=sample_std)
print(conf_intveral)
return conf_intveral
# T分布下的置信区间
def ttest_conf2(data, confidence=0.95):
# 先求一下 bins ,以便画图用得上.
data2 = pd.DataFrame(data, columns=['age'])
IQR = data2.quantile(0.75) - data2.quantile(0.25)
bin_size =2 * IQR /len(data)**(1.0/3)
# 画个源数据图表.以便对源数据的一个直观了解
plt.rcParams['font.sans-serif']=['SimHei']
n = plt.hist(data,bins =round(bin_size),rwidth=0.9)
plt.vlines(data2.mean(),0,max(n[1])+1,colors ="r", linestyles ="dashed",label="平均值%.2f" % np.mean(data))
plt.title('源数据'+str(len(data))+'个样本分布 直方图')
plt.ylabel('频数')
plt.legend()
plt.show()
# 真正开始计算
sample_mean = np.mean(data2)
sample_std = np.std(data2)
sample_size =len(data2)
alpha =1 -0.95
t_score = scipy.stats.t.isf(alpha /2, df = (sample_size-1) )
ME = t_score * sample_std / np.sqrt(sample_size)
lower_limit = sample_mean - ME
upper_limit = sample_mean + ME
print(str(confidence*100)+'%% Confidence Interval: ( %.2f, %.2f)' % (lower_limit, upper_limit))
return lower_limit, upper_limit
if __name__ =="__main__":
try:#重复抽取数据1000次
scale_means = []
for _in range(1000):
scale_sample = data.sample(100, replace=True)
mean = scale_sample.mean()
scale_means.append(mean)
norm_conf(scale_means)#正态分布
ttest_conf(scale_means)#t分布
ttest_conf2(scale_means)#t分布
import seabornas sns
from matplotlibimport pyplotas plt
sns.set_palette("hls")# 设置所有图的颜色，使用hls色彩空间
sns.distplot(scale_means, color="r", bins=10, kde=True)
plt.title('Age')
plt.xlim(25, 35)
plt.grid(True)
plt.show()
except:
traceback.print_exc()
截图如下:
