2017-07-17 5 views
1

Ist es möglich, die Summe der Gewichte in jedem Fach eines Histogramms in numpy (oder scipy) zu erhalten? Ich möchte den Fehler auf jeder Behälterhöhe in meinen Histogrammen haben. Für nicht gewogene Daten sollte der statistische Fehler auf jeder Behälterhöhe der Wert sqrt (N) sein, wobei N die Behälterhöhe ist. Für gewichtete Daten muss ich jedoch die Summe der Gewichte im Quadrat berechnen. numpy.histogram kann dies nicht tun, aber gibt es eine andere Funktionalität in numpy oder scipy, die ein Array (z. B. die Gewichtungsmatrix) basierend auf einem anderen Array (z. B. das Array von Werten, die ich histogramming) bin? Ich habe die Dokumentation gelesen, aber nichts gefunden.numpy.histogram: Gewichtete Summe der Gewichte in jedem Fach abfragen

+2

Ich würde beginnen mit [numpy.digitize] (https://docs.scipy.org/doc/numpy/reference/generated/numpy.digitize.html#numpy.digitize) – FTP

+0

Ich verstehe es nicht. Könnten Sie das mathematischer ausdrücken? – obachtos

+0

@obachtos Sagen wir, ich habe ein Array 'x = [2,9,4,8]' und ein Array von Gewichten 'w = [0,1,0,2,0,3,0,4] .' Ich erstelle ein Histogramm mit 2 Bins mit 'numpy.histogram (x, gewichte = w, bins = [0,5,10])'. In der 0-ten Ablage bekomme ich die 2 und die 4, aber eine Gesamtkorbhöhe von 0,1 + 0,3 = 0,4 aufgrund der Gewichte. Im ersten bin ich die 9 und die 8 mit einer Behälterhöhe von 0,2 + 0,4 = 0,6. Ich möchte auch die Summe der Gewichte für jeden Bin_ quadriert bekommen. Das wäre .1^2 + .3^2 für den 0. Fach. der statistische Fehler auf dieser Behälterhöhe wäre sqrt (Summe (.1^2 + .3^2)) = 0.316 ... nicht sqrt (Behälterhöhe) wie bei ungewichteten Daten. – ddavis

Antwort

2

Wie Alex vorgeschlagen hat, ist numpy.digitize was Sie wollen. Diese Funktion gibt an, zu welchen Bins die Einträge Ihres Arrays x gehören. Sie können dann diese Informationen verwenden, um die richtigen Elemente von w zuzugreifen:

x = np.array([2,9,4,8]) 
w = np.array([0.1,0.2,0.3,0.4]) 

bins = np.digitize(x, [0,5,10]) 

# access elements for first bin 
first_bin_ws = w[np.where(bins==1)[0]] 

# error of fist bin 
error = np.sqrt(np.sum(first_bin_ws**2.)) 

Die letzte Zeile berechnet dann den Fehler für den ersten Behälter. Daran, dass np.digitize beginnt bei 1.

+0

perfekt - das Beispiel mit 'np.where' ist, was ich sehen musste, um es zu klicken. – ddavis

1

zu zählen Wenn ich eine Ergänzung hinzufügen können @ obachtos Antwort, habe ich es in eine Funktion erweitert, die dann für das gesamte Histogramm zeigt:

def hist_bin_uncertainty(data, weights, bin_edges): 
    """ 
    The statistical uncertainity per bin of the binned data. 
    If there are weights then the uncertainity will be the root of the 
    sum of the weights squared. 
    If there are no weights (weights = 1) this reduces to the root of 
    the number of events. 

    Args: 
     data: `array`, the data being histogrammed. 
     weights: `array`, the associated weights of the `data`. 
     bin_edges: `array`, the edges of the bins of the histogram. 

    Returns: 
     bin_uncertainties: `array`, the statistical uncertainity on the bins. 

    Example: 
    >>> x = np.array([2,9,4,8]) 
    >>> w = np.array([0.1,0.2,0.3,0.4]) 
    >>> edges = [0,5,10] 
    >>> hist_bin_uncertainty(x, w, edges) 
    array([ 0.31622777, 0.4472136 ]) 
    >>> hist_bin_uncertainty(x, None, edges) 
    array([ 1.41421356, 1.41421356]) 
    >>> hist_bin_uncertainty(x, np.ones(len(x)), edges) 
    array([ 1.41421356, 1.41421356]) 
    """ 
    import numpy as np 
    # Bound the data and weights to be within the bin edges 
    in_range_index = [idx for idx in range(len(data)) 
         if data[idx] > min(bin_edges) and data[idx] < max(bin_edges)] 
    in_range_data = np.asarray([data[idx] for idx in in_range_index]) 

    if weights is None or np.array_equal(weights, np.ones(len(weights))): 
     # Default to weights of 1 and thus uncertainty = sqrt(N) 
     in_range_weights = np.ones(len(in_range_data)) 
    else: 
     in_range_weights = np.asarray([weights[idx] for idx in in_range_index]) 

    # Bin the weights with the same binning as the data 
    bin_index = np.digitize(in_range_data, bin_edges) 
    # N.B.: range(1, bin_edges.size) is used instead of set(bin_index) as if 
    # there is a gap in the data such that a bin is skipped no index would appear 
    # for it in the set 
    binned_weights = np.asarray(
     [in_range_weights[np.where(bin_index == idx)[0]] for idx in range(1, len(bin_edges))]) 
    bin_uncertainties = np.asarray(
     [np.sqrt(np.sum(np.square(w))) for w in binned_weights]) 
    return bin_uncertainties