python 计算数据中值的置信区间

2019-06-03  本文已影响0人  YuAllon

最近给导师报告处理WISE数据的进展时,对于处理WISE光变导师希望我对每一个观测区间的星等数值给中值的置信度。

计算数据的中值95%的置信区间

假设有一组数据如下:

data = [-0.1, -2.4, -0.1, -0.7, -1.4, -0.9, -3.2, -0.2, -0.3, -0.6, -3.2, -5.5]

求中值的置信区间与求数据的均值置信区间的方法是类似的,在这里我不会讨论详细的数学原理,而是直接给出置信度为95%时,对应的数据的上限值与下限值。(有关python实现的重点在于代码)
下限:0.5n-0.98\sqrt{n}
上限:1+0.5n+0.98\sqrt{n}
上面公式给出的是理论值,具体应用到数据上时,要对得到的下限上限取整下限值向上取整上限值向下取整

一般情况:
Lower lim = n \over 2 - \sqrt{n} \over 2 * N_{1 - {\alpha \over 2}}
Upper lim = 1 + n \over 2 + \sqrt{n} \over 2 * N_{1 - {\alpha \over 2}}
注:n是数据个数,一般要求数据点个数n>=6。当n<6时是没有中值的置信区间。

计算上述data的中值95%置信区间

首先要将原数据从小到大排列:

sorted(data)
data1 = [-5.5,-3.2,-3.2,-2.4,-1.4,-0.9,-0.7,-0.6,-0.3,-0.2,-0.1,0.1]

下限值:0.5* 12 - 0.98* \sqrt{12} = 2.6
上限值:1+0.5* 12 + 0.98* \sqrt{12} = 10.4
则95%置信区间对应的数值是第3个数据第10个数据,即(-3.2,-0.2)

python实现

#求中值median的置信区间(confidence interval),95%的CI
#对于中值的置信区间CI,下限lower limit向上取整,upper limit向下取整
#要注意python中是从0开始计数的,根据上述就很好理解return语句的含义了。
#其实可以吧math.ceil() - 1 用math.floor()代替
import math
import numpy as np
import matplotlib.pyplot as plt

#自定义的median_ci函数是给出某一数据95%置信区间的上限和下限对应的值
def median_ci(data,confidence=0.95):
    data1 = sorted(data)
    n = len(data1)
    ll = 0.5*n - 0.98*math.sqrt(n)
    ul = 1 + 0.5*n + 0.98*math.sqrt(n)
    l = data1[math.ceil(ll)-1]
    u = data1[math.floor(ul) - 1]
    return (l,u)
#在自定义的函数里面已经将数据从小到大排序了,所以调用函数时用的是数据data
l,u = median_ci(data)
print(l,u)
-3.2,-0.2

参考:Confidence Intervals for a MedianDocumentation for
Confidence Interval of median or percentiles for a given sample size
.

上一篇下一篇

猜你喜欢

热点阅读