(Python) Szacowanie parametrów regresji przedziały ufności z scikits bootstrap

Question

Właśnie zaczął wypróbować piękny pakiet ładowania dostępne poprzez scikits: https://github.com/cgevans/scikits-bootstrap

ale ja napotkał problem podczas próby oszacowania przedziały ufności dla współczynnik korelacji z regresji liniowej. Zwarte przedziały ufności leżą całkowicie poza zakresem oryginalnych statystyk.

tutaj kod:

import numpy as np 
from scipy import stats 
import bootstrap as boot 

np.random.seed(0) 
x  = np.arange(10) 
y  = 10 + 1.5*x + 2*np.random.randn(10) 
r0 = stats.linregress(x, y)[2] 

def my_function(y): 
    return stats.linregress(x, y)[2] 

ci = boot.ci(y, statfunction=my_function, alpha=0.05, n_samples=1000, method='pi') 

Daje to w wyniku [CI = -0.605, 0,644], ale przy parametrem jest R0 = 0,894.

Próbowałem to w R i wydaje się, że działa dobrze tam: ci okrakiem r0 zgodnie z oczekiwaniami.

Proszę pomóc!

Answer 1

Czy możesz podać swój kod R? Byłbym zainteresowany wiedząc, jak to omówiono w R.

Problem polega na tym, że jesteś tylko przechodząc y do boot.ci, ale za każdym razem przebiega My_function, wykorzystuje on cały, oryginalny x (zauważ brak x wejścia do funkcji my_function). Bootstrapping stosuje funkcję statystyczną do resamplingowania danych, więc jeśli zastosujesz funkcję statystyczną używając oryginalnego x i próbki y, uzyskasz bezsensowny wynik. Właśnie dlatego metoda BCA w ogóle nie działa: nie może zastosować twojej statystyki do próbek scyzoryków, które nie mają takiej samej liczby elementów.

Próbki są pobierane wzdłuż osi 0 (wiersze), więc jeśli chcesz przekazać wiele tablic 1D do swojej funkcji statystycznej, możesz użyć wielu kolumn: xy = vstack((x,y)).T będzie działać, a następnie użyć funkcji statystycznej, która pobiera dane z tych kolumn:

def my_function(xysample): 
    return stats.linregress(xysample[:,0], xysample[:,1])[2] 

Ewentualnie, jeśli chciał uniknąć brudząc z danymi w ogóle, można zdefiniować funkcję, która działa na indeksach, a potem po prostu przejść indeksy boot.ci:

def my_function2(i): 
    return stats.linregress(x[i], y[i])[2] 

boot.ci(np.arange(len(x)), statfunction=my_function2, alpha=0.05, n_samples=1000, method='pi') 

Należy pamiętać, że w każdym z tych przypadków, BCA działa, więc równie dobrze możesz użyć metody = 'bca', chyba że naprawdę chcesz korzystać z interwałów procentowych; BCA jest prawie zawsze lepszy.

Zdaję sobie sprawę, że obie te metody są mniej niż idealne. Szczerze mówiąc, nigdy nie musiałem przekazywać wielu tablic takich jak ta do mojej funkcji statycznej, a większość ludzi prawdopodobnie używa mean jako ich funkcji statfunkcji. Myślę, że najlepszym pomysłem może być tutaj umożliwienie przekazywania list o równych rozmiarach [0], np. boot.ci([x,y],...), a następnie wypróbowanie wszystkich tych tablic w tym samym czasie i przekazanie ich wszystkich do statfunkcji jako oddzielnych argumentów. W takim przypadku możesz po prostu mieć my_function(x,y). Zobaczę, czy mogę to zrobić, ale jeśli możesz pokazać mi swój kod R, to byłoby wspaniale, ponieważ chciałbym zobaczyć, czy jest lepszy sposób radzenia sobie z tym.

---

Aktualizacja:

W najnowszej wersji scikits.bootstrap (v0.3.1), krotka tablic mogą być zrealizowane, a próbki z nich zostanie przekazana jako osobne argumenty statfunction. Dodatkowo, funkcja statfunction może dostarczyć danych wyjściowych tablicy, a przedziały ufności zostaną obliczone dla każdego punktu na wyjściu.Jest to teraz bardzo łatwe do zrobienia. Następujące daje przedziały ufności dla każdego wyjścia linregress:

cis = boot.ci((x,y), statfunction=stats.linregress) 

cis[:,2] w tym przypadku będzie pożądany przedział ufności.