Рассчитать FDR вручную в Python - PullRequest
Новая
       72

Рассчитать FDR вручную в Python

0 голосов
/ 21 декабря 2018

Я пытаюсь вручную вычислить частоту ложных обнаружений (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 <- sort(p)
fdrtest <- NULL
for (i in 1:1000)
    fdrtest <- c(fdrtest, p[i] > 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

мой тип IIуровень ошибок намного выше (ближе к 60%) по сравнению с примером ближе к 40%.Ниже показано, как я вычисляю ошибки типа 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)

Справочный документ

...