Я пытаюсь написать функцию, которая возвращает значение для оценки пробега шкалы набора данных. Может кто-нибудь взглянуть?

Question

Если вам не нужен контекст и вы просто хотите посмотреть на код, пропустите два абзаца.

Я работаю над проектом, где я измеряю дисперсии скорости скоплений галактик, учитывая данные красного смещения галактики, собранные с телескопа нашей группы на южном полюсе (SPT). Для кластеров, которые имеют очень мало галактик-членов, типичных оценок, таких как стандартное отклонение или даже более сложная дисперсия в весовом диапазоне, недостаточно.

Wainer & Thissen (1976) предлагает интересную оценку масштаба, которая учитывает разрывы между данными по своей мере дисперсии, которые они показывают, что они очень эффективны для наборов данных небольших размеров (93% при n = 10)

Вот определение:

, где п является размер выборки, а

, промежутки между точками данных

, набор примерно гауссовских весов и

Его довольно просто, но я пытаюсь реализовать в Python и продолжать получать неправильно ответ. Ответы слишком малы. Я не вижу ничего плохого в этом, хотя, я проверял почти все. Кто-нибудь видит тонкую синтаксическую ошибку, которую я пропускаю или что-то еще? Вот эта функция:

def gDispersion(v): 
    # Returns the gapper velocity dispersion of a cluster (Sigma_G), given galaxy proper velocity data, 
    # v is an array of velocity values. 

    try: 
     #allocate array for Gaussian weights 
     w = [0] * (len(v)-1) 
     g = [0] * (len(v)-1) 
     n = len(v) 
    except TypeError: 
     # ensure input is valid 
     print('Array or array-like object expected; got {}'.format(type(v))) 
     return 

    # find gaussian weights 
    for i in range(len(v) - 1): 
     g[i] = v[i+1] - v[i] 
     w[i] = i * (n - i) 

    sigG = (np.sqrt(np.pi))/(n*(n-1)) * sum([wi*gi for wi,gi in zip(w,g)]) 
    return sigG 

Answer 1

Ваши весы неправильные. В цикле i принимает значения 0, 1, ..., n-2, так что

w[0] = 0, w[1] = 1*(n-1), ..., w[n-2] = (n-2)*2 

Цикл также ненужным. Вот Векторизованная реализация:

def gDispersion(v): 
    """ Returns the gapper velocity dispersion of a cluster (Sigma_G) 

    v is an array of galaxy velocity values. 
    """ 
    n = len(v) 
    w = np.arange(1, n) * np.arange(n-1, 0, -1) 
    g = np.diff(v) 
    sigG = (np.sqrt(np.pi))/(n*(n-1)) * np.dot(w, g) 
    return sigG 

Кстати, использовать docstrings вместо комментариев к документу своих функций.