python手工制作_用python手工计算fdr

import scipy.stats as sc_stats

import numpy as np

np.random.seed(471)

a = np.random.normal(0, 1, 900)

b = np.random.normal(3, 1, 100)

x = np.concatenate([a, b])

p = sc_stats.norm.pdf(x)

根据样本计算FDR的方法如下:

对于fdr,我们要考虑有序p值。我们来看看第k个有序p值是否大于k*0.05/1000

psort

fdrtest

for (i in 1:1000)

fdrtest match(p[i],psort) * .05/1000)

我已经使用一个详细的函数实现了如下:

def get_index(s, item):

for i, si in enumerate(s): if si >= item: return i

return i

psort = np.sort(p)

fdrtest = np.zeros(len(p))

wheres = np.zeros(len(p))

for i, pi in enumerate(p):

#print(i, np.where(psort==pi)[0][0], pi, psort[np.where(psort==pi)[0]])

wheres[i] = get_index(psort, pi)

fdrtest[i] = True if pi > (wheres[i] * 0.05/1000) else False

与接近40%的示例相比,我的II型错误率要高得多(接近60%)下面是我如何计算类型I和类型II错误。

tmp = fdrtest

a_test = pd.DataFrame({'vals':

fdrtest[:900]}).groupby('vals').size().reset_index()

display(a_test)

b_test = pd.DataFrame({'vals':

fdrtest[900:]}).groupby('vals').size().reset_index()

display(b_test)

a_reject_ct = a_test[a_test['vals']==False].values[0][1]

a_samples = a.shape[0]

b_samples = b.shape[0]

b_fail_rej_ct = b_test[b_test['vals']==True].values[0][1]

type_1_err = np.round(a_reject_ct / a_samples,4)

type_2_err = np.round( b_fail_rej_ct/ b_samples,4)