Ajustement d'un histogramme avec python

Question

J'ai un histogramme

H = hist(mon_data, bins = mon_bin, histtype = 'step', color = 'r')

Je peux voir que la forme est presque gaussienne mais j'aimerais ajuster cet histogramme avec une fonction gaussienne et afficher les valeurs de la moyenne et de sigma que j'obtiens. Pouvez-vous m'aider?

Answer 1

Voici un exemple fonctionnant sur py2.6 et py3.2 :

from scipy.stats import norm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt

# lire les données à partir d'un fichier texte. Un nombre par ligne
arch = "test/Log(2)_ACRatio.txt"
donnees = []
for item in open(arch,'r'):
    item = item.strip()
    if item != '':
        try:
            donnees.append(float(item))
        except ValueError:
            pass

# meilleur ajustement des données
(mu, sigma) = norm.fit(donnees)

# l'histogramme des données
n, bins, patches = plt.hist(donnees, 60, normed=1, facecolor='green', alpha=0.75)

# ajouter une ligne de 'meilleur ajustement'
y = mlab.normpdf( bins, mu, sigma)
l = plt.plot(bins, y, 'r--', linewidth=2)

# tracer
plt.xlabel('Intelligence')
plt.ylabel('Probabilité')
plt.title(r'$\mathrm{Histogramme\ de\ QI:}\ \mu=%.3f,\ \sigma=%.3f$' %(mu, sigma))
plt.grid(True)

plt.show()

Answer 2

Voici un exemple qui utilise scipy.optimize pour ajuster des fonctions non linéaires comme une Gaussienne, même lorsque les données sont dans un histogramme qui n'est pas bien étalonné, de sorte qu'une simple estimation de la moyenne échouerait. Une constante de décalage causerait également l'échec de statistiques normales simples (il suffit de supprimer p[3] et c[3] pour des données gaussiennes simples).

from pylab import *
from numpy import loadtxt
from scipy.optimize import leastsq

fitfunc  = lambda p, x: p[0]*exp(-0.5*((x-p[1])/p[2])**2)+p[3]
errfunc  = lambda p, x, y: (y - fitfunc(p, x))

filename = "gaussdata.csv"
data     = loadtxt(filename,skiprows=1,delimiter=',')
xdata    = data[:,0]
ydata    = data[:,1]

init  = [1.0, 0.5, 0.5, 0.5]

out   = leastsq( errfunc, init, args=(xdata, ydata))
c = out[0]

print "A exp[-0.5((x-mu)/sigma)^2] + k "
print "Coefficients Parentaux:"
print "1.000, 0.200, 0.300, 0.625"
print "Coefficients Ajustés:"
print c[0],c[1],abs(c[2]),c[3]

plot(xdata, fitfunc(c, xdata))
plot(xdata, ydata)

title(r'$A = %.3f\  \mu = %.3f\  \sigma = %.3f\ k = %.3f $' %(c[0],c[1],abs(c[2]),c[3]));

show()

Résultat:

A exp[-0.5((x-mu)/sigma)^2] + k 
Coefficients Parentaux:
1.000, 0.200, 0.300, 0.625
Coefficients Ajustés:
0.961231625289 0.197254597618 0.293989275502 0.65370344131

Answer 3

À partir de Python 3.8, la bibliothèque standard fournit l'objet NormalDist en tant que partie du module statistics.

L'objet NormalDist peut être construit à partir d'un ensemble de données avec la méthode NormalDist.from_samples et donne accès à sa moyenne (NormalDist.mean) et à son écart type (NormalDist.stdev):

from statistics import NormalDist

# data = [0.7237248252340628, 0.6402731706462489, -1.0616113628912391, -1.7796451823371144, -0.1475852030122049, 0.5617952240065559, -0.6371760932160501, -0.7257277223562687, 1.699633029946764, 0.2155375969350495, -0.33371076371293323, 0.1905125348631894, -0.8175477853425216, -1.7549449090704003, -0.512427115804309, 0.9720486316086447, 0.6248742504909869, 0.7450655841312533, -0.1451632129830228, -1.0252663611514108]
norm = NormalDist.from_samples(data)
# NormalDist(mu=-0.12836704320073597, sigma=0.9240861018557649)
norm.mean
# -0.12836704320073597
norm.stdev
# 0.9240861018557649

Answer 4

Voici une autre solution utilisant uniquement les packages matplotlib.pyplot et numpy. Cela fonctionne uniquement pour l'ajustement gaussien. Il est basé sur l'estimation du maximum de vraisemblance et a déjà été mentionné dans ce sujet. Voici le code correspondant :

# Version de Python : 2.7.9
from __future__ import division
import numpy as np
from matplotlib import pyplot as plt

# Pour l'explication, je simule les données :
N=1000
data = np.random.randn(N)
# Mais en réalité, vous liriez les données à partir d'un fichier, par exemple avec :
#data = np.loadtxt("data.txt")

# La moyenne empirique et la variance sont calculées
avg = np.mean(data)
var = np.var(data)
# À partir de là, nous connaissons la forme de la gaussienne ajustée.
pdf_x = np.linspace(np.min(data),np.max(data),100)
pdf_y = 1.0/np.sqrt(2*np.pi*var)*np.exp(-0.5*(pdf_x-avg)**2/var)

# Ensuite, nous traçons :
plt.figure()
plt.hist(data,30,normed=True)
plt.plot(pdf_x,pdf_y,'k--')
plt.legend(("Ajustement","Données"),"best")
plt.show()

et voici la sortie.

Answer 5

J'étais un peu perplexe que norm.fit ne semblait fonctionner qu'avec la liste étendue des valeurs échantillonnées. J'ai essayé de lui donner deux listes de nombres, ou des listes de tuples, mais il semblait tout aplatir et traiter l'entrée comme des échantillons individuels. Comme j'ai déjà un histogramme basé sur des millions d'échantillons, je ne voulais pas l'étendre si je n'y étais pas obligé. Heureusement, la distribution normale est facile à calculer, donc...

# histogramme est [(val,count)]
from math import sqrt

def normfit(hist):
    n,s,ss = univar(hist)
    mu = s/n
    var = ss/n-mu*mu
    return (mu, sqrt(var))

def univar(hist):
    n = 0
    s = 0
    ss = 0
    for v,c in hist:
        n += c
        s += c*v
        ss += c*v*v
    return n, s, ss

Je suis sûr que cela doit être fourni par les bibliothèques, mais comme je ne l'ai trouvé nulle part, je le poste ici à la place. N'hésitez pas à indiquer la meilleure façon de le faire et à me déclasser :-)