Contabilizzazione degli errori durante la creazione di un istogramma

Question

Ho un set di osservazioni N distribuite come punti (x[i], y[i]), i=0..N in uno spazio 2D. Ogni punto ha errori associati in entrambe le coordinate (e_x[i], e_y[i], i=0..N) e anche un peso collegato ad esso (w[i], i=0..N).

mi piacerebbe per generare un istogramma 2D di questi N punti, che rappresentano non solo per i pesi, ma anche per gli errori, che causerebbe ogni punto di essere diffondere possibilmente tra molti bidoni se i valori di errore sono grandi abbastanza (assumendo uno standard Gaussian distribution per gli errori, anche se potrebbero essere prese in considerazione altre distribuzioni).

Vedo che numpy.histogram2d ha un parametro weights che viene preso in considerazione. Il problema sarebbe come tenere conto degli errori in ciascuno dei punti osservati N.

Esiste una funzione che consenta di eseguire questa operazione? Sono aperto a qualsiasi cosa in numpy e scipy.

Answer 1

Basandosi sul commento dell'utente 1415946, è possibile assumere che ogni punto rappresenti uno bi-variate normal distribution con le matrici di covarianza fornite da [[e_x[i]**2,0][0,e_y[i]**2]]. Tuttavia, la distribuzione risultante non è una distribuzione normale - vedrai, dopo aver eseguito l'esempio, come l'istogramma non assomigli affatto a un gaussiano, ma piuttosto un gruppo di essi.

Per creare un istogramma da questo insieme di distribuzioni, un modo che vedo è generare campioni casuali su ciascun punto utilizzando numpy.random.multivariate_normal. Vedi il codice di esempio qui sotto con alcuni dati artificiali.

import numpy as np 
from mpl_toolkits.mplot3d import Axes3D 
import matplotlib.pyplot as plt 


# This is a function I like to use for plotting histograms 
def plotHistogram3d(hist, xedges, yedges): 
    fig = plt.figure() 
    ax = fig.add_subplot(111, projection='3d') 
    hist = hist.transpose() 
    # Transposing is done so that bar3d x and y match hist shape correctly 
    dx = np.mean(np.diff(xedges)) 
    dy = np.mean(np.diff(yedges)) 

    # Computing the number of elements 
    elements = (len(xedges) - 1) * (len(yedges) - 1) 
    # Generating mesh grids. 
    xpos, ypos = np.meshgrid(xedges[:-1]+dx/2.0, yedges[:-1]+dy/2.0) 

    # Vectorizing matrices 
    xpos = xpos.flatten() 
    ypos = ypos.flatten() 
    zpos = np.zeros(elements) 
    dx = dx * np.ones_like(zpos) * 0.5 # 0.5 factor to give room between bars. 
# Use 1.0 if you want all bars 'glued' to each other 
    dy = dy * np.ones_like(zpos) * 0.5 
    dz = hist.flatten() 

    ax.bar3d(xpos, ypos, zpos, dx, dy, dz, color='b') 
    ax.set_xlabel('x') 
    ax.set_ylabel('y') 
    ax.set_zlabel('Count') 
    return 

""" 
INPUT DATA 
""" 
#     x y ex ey w 
data = np.array([[1, 2, 1, 1, 1], 
       [3, 0, 1, 1, 2], 
       [0, 1, 2, 1, 5], 
       [7, 7, 1, 3, 1]]) 

""" 
Generate samples 
""" 
# Sample size (100 samples will be generated for each data point) 
SAMPLE_SIZE = 100 
# I want to fill in a table with columns [x, y, w]. Each data point generates SAMPLE_SIZE 
# samples, so we have SAMPLE_SIZE * (number of data points) generated points 
points = np.zeros((SAMPLE_SIZE * data.shape[0], 3)) # Initializing this matrix 

for i, element in enumerate(data): # For each row in the data set 
    meanVector = element[:2] 
    covarianceMatrix = np.diag(element[2:4]**2) # Diagonal matrix with elements equal to error^2 
    # For columns 0 and 1, add generated x and y samples 
    points[SAMPLE_SIZE*i:SAMPLE_SIZE*(i+1), :2] = \ 
     np.random.multivariate_normal(meanVector, covarianceMatrix, SAMPLE_SIZE) 
    # For column 2, simply copy original weight 
    points[SAMPLE_SIZE*i:SAMPLE_SIZE*(i+1), 2] = element[4] # weights 

hist, xedges, yedges = np.histogram2d(points[:, 0], points[:, 1], weights=points[:, 2]) 
plotHistogram3d(hist, xedges, yedges) 
plt.show() 

Risultati tracciate di seguito: