TUTORIAL PROFILE PLOT¶

In questo tutorial viene mostrato come costruire in modo rapido ed indolore un profile plot usando il pitone.

Lo stimolo è nato da una problematica avuta in laboratorio. Il punto di partenza è l'esempio di Valerio da cui prenderò anche ispirazione per i dati da usare. Il motivo è che tutte le funzioni che ho visto utilizzano dei cicli for che, oltre ad essere il male su Python, non sono veramente necessari ed anzi rischiano di rallentare il codice e di renderlo meno comprensibile

In [1]:
import numpy as np
from matplotlib import pyplot as plt

Genero e plotto dei dati casuali. Utilizzo un differente binnaggio per poter seguire al meglio quello che succede, e non trovarmi a dover riflettere se lavorare sulla matrice hh[0] o sulla sua trasposta. Lungo la prima dimensione (ovvero muovendosi lungo le righe) viene riportata la suddivisione in x, che nel plot viene mostrata in orizzontale

In [2]:
N = 10000
x = np.random.rand(N)*2*np.pi # flat distribution between 0 and 2pi
y = np.sin(x) + np.random.normal(0, 0.1, N) +  np.random.normal(0, 0.3, N)*x/5 # sin wave + some spread which increases with x 

fig, ax = plt.subplots()
mat, binsx, binsy, hh = ax.hist2d(x,y, (50, 60), cmap='coolwarm')
fig.colorbar(hh, ax = ax)
plt.show()
In [3]:
print(mat.shape, binsx.shape, binsy.shape, )

bincx = binsx[:-1] + (binsx[1]-binsx[0])/2
bincy = binsy[:-1] + (binsy[1]-binsy[0])/2

print(mat.shape, bincx.shape, bincy.shape, )

(50, 60) (51,) (61,)
(50, 60) (50,) (60,)

Per costruire il profile plot devo prendere, per ciascuna colonna del plot sopra, valor medio e deviazione standard. Avendo scelto 50 valori di x, mi aspetto di avere 50 punti. Quando uso la notazione tile in questo modo, la dimensione del vettore originale diventa la seconda

In [4]:
maty = np.tile(bincy, (bincx.shape[0],1))
print(maty.shape)

(50, 60)
  • Mi creo la media come media pesata su ciascuna colonna, essendo i pesi i conteggi
  • Mi creo una matrice che mi replica per ciascuna riga il vettore delle medie (in pratica tutte le colonne sono uguali)
  • Calcolo la standard deviazione (cfr. questo link)
In [5]:
prodotto = maty * mat
print(prodotto.shape)

# media pesata ed errore sulla media pesata
myMean = np.nansum(prodotto, axis = 1) / np.nansum(mat, axis=1)
myMeanMat = np.tile(myMean, (bincy.shape[0],1)).T

myStd = np.sqrt( np.sum( mat*((maty-myMeanMat)**2), axis = 1 ) / (np.sum(~np.isnan(mat), axis = 1) - 1))


print(myMean.shape, myStd.shape)
print(myMeanMat.shape)

(50, 60)
(50,) (50,)
(50, 60)

Disegno il profile plot

In [6]:
fig, ax = plt.subplots()

ax.errorbar(bincx, myMean, yerr = myStd, marker = "*", c="tab:green", ls = "")

plt.show()

Mettendo tutto assieme¶

In [7]:
# Genero la distrubuzione dei dati
N = 10000
x = np.random.rand(N)*2*np.pi # flat distribution between 0 and 2pi
y = np.sin(x) + np.random.normal(0, 0.1, N) +  np.random.normal(0, 0.3, N)*x/5 # sin wave + some spread which increases with x 



fig, ax = plt.subplots()

# Istogrammo i dati
mat, binsx, binsy, hh = ax.hist2d(x,y, (50, 60), cmap='coolwarm')

# Centri dei bins
bincx = binsx[:-1] + (binsx[1]-binsx[0])/2
bincy = binsy[:-1] + (binsy[1]-binsy[0])/2


# Ottengo una matrice con i bin in y
maty = np.tile(bincy, (bincx.shape[0],1))

# media pesata ed errore sulla media pesata
myMean = np.nansum(maty * mat, axis = 1) / np.nansum(mat, axis=1)
myMeanMat = np.tile(myMean, (bincy.shape[0],1)).T

myStd = np.sqrt( np.sum( mat*((maty-myMeanMat)**2), axis = 1 ) / (np.sum(~np.isnan(mat), axis = 1) - 1))


# Plotto il profile
ax.errorbar(bincx, myMean, yerr = myStd, marker = ".", c = "k", ls = "")




fig.colorbar(hh, ax = ax)
plt.show()

Definisco due funzioni comode¶

  • (Creare istogramma e profile)
  • Creare profile a partire da un istogramma

// TODO: Implementare documentazione seria e varie possibilità per argomenti

In [8]:
def istoAndProfile(x, y, **kwargs):
    pass
In [9]:
def profileFromIsto(isto):
    mat = isto[0]
    binsx = isto[1]
    binsy = isto[2]
    hh = isto[3]
    
    
    # Centri dei bins
    bincx = binsx[:-1] + (binsx[1]-binsx[0])/2
    bincy = binsy[:-1] + (binsy[1]-binsy[0])/2


    # Ottengo una matrice con i bin in y
    maty = np.tile(bincy, (bincx.shape[0],1))

    # media pesata ed errore sulla media pesata
    myMean = np.nansum(maty * mat, axis = 1) / np.nansum(mat, axis=1)
    myMeanMat = np.tile(myMean, (bincy.shape[0],1)).T

    myStd = np.sqrt( np.sum( mat*((maty-myMeanMat)**2), axis = 1 ) / (np.sum(~np.isnan(mat), axis = 1) - 1))
    
    return (bincx, myMean, myStd)
In [10]:
fig, ax = plt.subplots()
myIsto = ax.hist2d(x,y, (50, 60), cmap='coolwarm')

# Chiamo la funzione passando l'oggetto istogramma
# Lei mi ritorna x, y ed errori
bincx, myMean, myStd = profileFromIsto(myIsto)
ax.errorbar(bincx, myMean, yerr = myStd, marker = ".", c = "k", ls = "")


plt.show()
In [ ]: