J'essaye d'écrire une fonction qui renvoie la valeur pour une estimation de l'échelle d'un ensemble de données. Quelqu'un peut-il jeter un coup d'oeil?

Question

Si vous ne vous souciez pas du contexte et que vous voulez simplement regarder le code, sautez deux paragraphes.

Je travaille sur un projet où je mesure des dispersions de vitesse de grappes de galaxies, étant donné les données de redshift de la galaxie recueillies à partir du télescope de notre groupe au pôle sud (SPT). Pour les grappes qui ont très peu de galaxies membres, les estimateurs typiques tels que l'écart-type ou même une dispersion plus fine de biweight ne sont pas suffisants. Thissen (1976) propose un estimateur d'échelle intéressant qui considère les écarts entre les données dans sa mesure de la dispersion, qu'ils montrent très efficace pour les ensembles de données de petites tailles (93% à n = 10)

Voici la définition:

où n est la taille de l'échantillon, et

, les écarts entre les points de données

, un ensemble de poids approximativement gaussienne, et

Son assez simple, mais je suis en train de le mettre en œuvre en Python et continuer à obtenir le mauvais répondre. Des réponses qui sont trop petites. Je ne vois rien de mal à cela, j'ai presque tout inspecté. Est-ce que quelqu'un voit une erreur de syntaxe subtile qui me manque ou quelque chose? Voici la fonction:

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

Vos poids sont erronés. Dans la boucle i prend les valeurs 0, 1, ..., n-2, de sorte

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

La boucle est également inutile. Voici une mise en œuvre vectorisé:

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 

Par ailleurs, utiliser docstrings au lieu de commentaires pour documenter vos fonctions.