Graphique avec scipy ? matplotlib ?

Version imprimable

16/04/2013, 14h46
Isabella83

1 pièce(s) jointe(s)

Graphique avec scipy ? matplotlib ?
Bonjour à tous,
Je dois recréer le même type de graphique que l'image ci dessous.
Ici, les données représentent la distance relative d'une entité par rapport à un point. On voit que les entités sont la plupart du temps sur le point étudié.
En dessous, les petits carrés représentent chaque entité, c'est une autre façon de représenter la courbe.
J'ai plusieurs question :
1/ Comment puis je couper ma fenêtre graphique en 2, de façon à avoir au dessus, la courbe, et en dessous, les petits carrés représentant les entités ?
2/ Comment réaliser ces petits carrés ?

Pour les entités se trouvant sur le point par exemple, je pensais faire un fichier contenant les x et y que je souhaite plotter, et incrémenter le y de façon à ce que les points ne soit pas tous les un sur les autres quand j'ai plusieurs point ayant le même x
Donc, dans l'image, ci dessous, par exemple, le pic en position 0 a 174 entités donc faire :
Code:

1 2 3 4 5 6 7 8 0 0.5 0 1 . . . . 0 87
Comment faire en sorte que mes points soient représentés par des carrés ?

Merci d'avance pour votre aide,
Isa
16/04/2013, 15h55
Matthieu Brucher

1) pas la peine, affiche l'autre graphique en inversant y !
2) tu peux specifier le marqueur dans matplotlib
16/04/2013, 16h00
Isabella83

Ok ok je vais voir ça mais je n'y avais pas pensé :)

Salut.
Je donne un exemple d'une figure avec 2 zones graphiques :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 import numpy as np import matplotlib.pyplot as plt y = np.random.normal( loc=50, scale=5, size=1e3) plt.subplot(2,1,1) classes = np.linspace(0, 100, 20) plt.hist(y,bins=classes) plt.subplot(2,1,2) classes = np.arange(0,100,1) plt.hist(y,bins=classes) plt.show()
Je comprends les représentions comme étant deux histogrammes. La seconde ayant des classes de largeur 1.
Certes, si tu veux représenter les données avec des carrées, ça n'est pas possible avec un histogramme. Il faut le faire comme un ensemble de points, ta façon de faire marche mais me semble un peu compliquée.
Pas compris le code de couleurs sur les carrés.

1 pièce(s) jointe(s)

Je vois ça comme ça :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
import numpy as np
import matplotlib.pyplot as plt
 
y = np.random.normal( loc=50, scale=10, size=500)
 
plt.subplot(2,1,1)
classes = np.linspace(0, 100, 20)
plt.hist(y,bins=classes)
 
plt.subplot(2,1,2)
classes = np.arange(0,100,1)
h, bin_edges = np.histogram(y,classes)
for dist, eff in zip(classes,h):
    for e in np.arange(eff):
        plt.plot(dist,e+1,'sr')
 
plt.xlim(0,100)
plt.show()

1 pièce(s) jointe(s)

Ah oui effectivement c'est bien plus simple de plotter point par point !
J'ai un peu modifié :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import pylab
 
y = np.random.normal( loc=50, scale=10, size=500)
 
ax1=plt.subplot(2,1,1)
classes = np.linspace(0, 100, 20)
plt.hist(y,bins=classes)
 
 
#plt.subplot(2,1,2)
ax2 = plt.subplot(212, sharex=ax1)
classes = np.arange(0,100,1)
h, bin_edges = np.histogram(y,classes)
for dist, eff in zip(classes,h): # distance dans mon histogramme + effectif dans chaque classe
	for e in np.arange(eff): # np.arrange(3) = array([0, 1, 2])
		#print "%d,%d : %d" %(dist,eff,e)
		plt.plot(dist,-e,'sg') # pour qu'il n'y ait pas d'overlappe
 
plt.xlim(0,100)
ax1.xaxis.tick_top()
ax2.axes.get_xaxis().set_visible(False)
ax2.axes.get_yaxis().set_visible(False)
plt.setp( ax2.get_xticklabels(), visible=False)
plt.setp( ax2.get_yticklabels(), visible=False)
plt.show()

j'ai mis ce que ça donne en pièce jointe :)

Je mettrai le résultat final avec mes données !

17/04/2013, 12h47
Matthieu Brucher

Pour le deuxieme plot, j'aurai plutot utilise scatterplot, sans dout eplus rapide, mais ca doit revenir a peu pres a la meme chose.
17/04/2013, 16h20
Isabella83

2 pièce(s) jointe(s)
Oui effectivement, j'ai fait avec scatterplot, si jamais j'ai à changer les couleurs ce sera plus pratique !
J'ai donc le schéma 20m1 avec mes data. Mais j'aimerai que la courbe soit lisse, j'ai donc fait une interpolation, mais j'obtiens la courbe en rouge avec le code suivant :
Code:

1 2 3 4 tck= interpolate.splrep(enreg['distance'],enreg['effectif']) xnew= np.linspace(enreg['distance'][0],enreg['distance'][-1], 1000) ynew = interpolate.splev(xnew,tck,der=0) plt.plot(xnew, ynew, 'r')
Que me conseillez vous pour faire une courbe lisse ? (sans les pics en fait )
17/04/2013, 16h27
Matthieu Brucher

Tu peux passer a scatter plot la liste des couleurs que tu veux mettre ;)

L'interpolation est un probleme ouvert. La courbe rouge est valide :/ Peut-etre en prenant une interpolation cubique ou avec un ordre faible ?
17/04/2013, 16h41
Isabella83

Citation:

Envoyé par Matthieu Brucher

Peut-etre en prenant une interpolation cubique ou avec un ordre faible ?

Le schema représente une interpolation cubic, j'ai essayé avec un ordre faible (k=1 ou 2) mais cela me donne la même chose ....
17/04/2013, 16h54
Matthieu Brucher

Vu que ca semble etre des distributions de probabilites, il faut peut-etre passer par une estimation de gaussiennes et les dessiner ?
17/04/2013, 17h13
__dardanos__

Quelles informations souhaites tu extraire ?

Si c'est la position, l'amplitude et la largeur des pics, l'interpolation par des splines ne me semble pas la bonne méthode.
Il faut plutôt approcher chaque pic par une forme théorique (ou une gaussienne faute de mieux).

La méthode consiste à approcher un premier pic, puis à le soustraire au signal et recommencer avec le signal obtenu jusqu'à traiter l'ensemble des pics.
Il faut fixer un seuil pour éliminer le bruit de fond.
17/04/2013, 17h16
Isabella83

En fait ce que j'aimerai, c'est avoir une courbe plus arrondie au niveau de ses pics, pas avoir ces pic triangulaires ....
Je viens de trouver ce lien http://jpktd.blogspot.fr/2009/03/usi...l-density.html
Peut être que cela correspond plus à ce que je souhaite faire ...
17/04/2013, 17h28
Matthieu Brucher

Je pense, c'est a ce genre de choses que je pensais. Tu essaies dans ce cas d'estimer une dentiste de proba consistant en une somme de loies normales.
17/04/2013, 17h32
__dardanos__

+1 pour le lien.

Je pense que la meilleure méthode est d'approcher chaque pic par une gaussienne et de sommer l'ensemble des gaussiennes obtenues pour former une approximation du signal original.
Les pics auront une forme bien lisse, et leurs caractéristiques seront connues avec précision.
17/04/2013, 17h37
Isabella83
Citation:

Envoyé par __dardanos__

Je pense que la meilleure méthode est d'approcher chaque pic par une gaussienne

Donc pour ce faire, je créer une liste avec la valeur en abscisse de chaque pic (comme dans le lien où les pics sont à -3 et à 3)

Mais comment faire l'estimation ? Je ne comprends pas bien la formule dans le deuxième cas du lien
Code:

1 2 3 plt.plot(ind, alpha * stats.norm.pdf(ind, loc=mlow) + (1-alpha) * stats.norm.pdf(ind, loc=mhigh), color="r", label='normal mix')
Car j'ai plus que 2 pics dans mes datas

1 pièce(s) jointe(s)

J'ai fait comme dans l'exemple du lien :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#!/usr/bin/env python
# -*- coding:Utf-8 -*- 
 
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import pylab
from scipy import interpolate
from scipy import stats
 
 
data = np.dtype( [('distance',np.int)])
enreg = np.loadtxt('tmp2.txt', dtype=data, skiprows=True)
 
gkde=stats.gaussian_kde(enreg['distance'])
 
ind= np.linspace(enreg['distance'][0],enreg['distance'][-1], 1000)
 
kdepdf = gkde.evaluate(ind)
 
plt.figure()
# plot histgram of sample
plt.hist(enreg['distance'], bins=20, normed=1)
# plot data generating density
plt.plot(ind, stats.norm.pdf(ind), color="r", label='DGP normal')
# plot estimated density
plt.plot(ind, kdepdf, label='kde', color="g")
plt.title('Kernel Density Estimation')
plt.legend()
plt.show()

Et j'obtiens l'image que j'ai mise en pièce jointe, qu'en pensez vous ?

17/04/2013, 18h15
Matthieu Brucher

Ca me parait pas mal pour une mixture de deux gaussiennes. J'en aurai mis 3 a priori, mais je ne sais pas ou ca se regle ni si ca se regle en fait.

Pour la formule, il s'agit de la mixture (un poids de alpha pour l'une et 1-alpha pour la seconde) de deux gaussiennes.
17/04/2013, 18h22
Isabella83

1 pièce(s) jointe(s)

Comment faire lorsque j'ai plus que 2 points où j'ai un pic ?
Car j'ai essayé en mettant 2 valeurs et j'obtiens l'image en pièce jointe, cela me semble étrange ...
18/04/2013, 08h29
__dardanos__

Salut.

J'ai déjà réalisé un travail analogue lors du développement d'un logiciel aidant à calibrer des détecteurs à particules. Le signal électronique enregistré par le détecteur présente des pics correspondant à l'énergie des particules interagissant avec le matériau. L'objectif est de déterminer avec précision la position des pics dans un spectre (le calibrage est possible car les énergies sont connues).
Les pics avaient la plupart du temps l'allure d'une gaussienne, parfois d'une gaussienne avec une "trainée". Le bruit de fond du signal étaient à prendre en compte lors du traitement.

Les spectres étaient composés de plusieurs pics. Il était nécessaire de traiter les pics par ordre d'énergie moyenne croissante.
Le traitement d'un pic est le suivant. Il faut estimer la fonction d'ajustement (par exemple par une gaussienne). Donner des valeurs initiales aux paramètres de la fonction (dans la cas général, pour une gaussienne 4 paramètres : position, largeur, amplitude et bruit de fond) et définir le domaine d'ajustement (les deux bornes de l'intervalle). J'utilise la fonction de scipy.interpolate.curve_fit. La fonction est totalement déterminée à partir des 4 paramètres.
Ayant la fonction d'ajustement, ce pic est soustrait au spectre initial. Le pic suivant dans le nouveau spectre est ensuite traité. Et ainsi de suite, jusqu'à traiter l'ensemble des pics.
Au final, un spectre initial comportant n pics est approché par un ensemble de gaussiennes définies par 4 paramètres chacune (autrement dit 4n paramètres suffisent à modéliser le spectre).

Une interface graphique avait été développée pour sélectionner à l'aide d'un Tcanvas de Tkinter (on pourrait faire la même chose avec une figure de Matplotlib) le domaine à traiter (et donc sélectionner un pic) et aider à définir les paramètres d'ajustement.

La méthode peut être reprise pour ce problème.

La première difficulté porte sur l'ensemble de points à approcher. Je comprends qu'il s'agit d'un histogramme de distances (je ne sais pas ce que représentent ces distances). Les résultats seront différents suivant la largeur des classes retenue. Plus la largeur est faible, meilleure est la précision (exemple de largeur valant 1 dans les figures du post 8), mais alors les effectifs des classes deviennent faibles.
Ainsi des pics peuvent apparaitre, ont-ils pour autant un sens physique ? Par exemple, sur les figures du post 8, je comprends qu'il y a trois pics entre les abscisses 15 et 25 (2 pics s'ajoutant à un pic beaucoup plus large). Alors qu'avec une largeur de classes plus importante, il n'y aurait qu'un seul. Comment alors choisir la largeur ?

La mise en œuvre est plus simple que dans mon problème car il n'y a pas de bruit de fond. Je traiterais alors les pics par ordre d'amplitude décroissante.
Il ne me semble pas indispensable de développer une interface graphique. Il est possible de définir à priori l'ordre de traitement. Sur les figures du post 8, d'abord le pic à 5, puis le pic situé autour de 20 (qui est mal défini avec cette largeur, la statistique est trop faible), et enfin le pic à -50 (mal défini aussi).

Bon courage.

Bonjour,
Merci pour toutes ces explications.
J'ai regardé du côté de curve_fit.
J'ai fait un test avec mes données :

Code:

et ce code :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#!/usr/bin/env python
# -*- coding:Utf-8 -*- 
 
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import pylab
from scipy import interpolate
from scipy.optimize import curve_fit
 
data = np.dtype( [ ('distance',np.int), ('effectif',np.int) ] )
enreg = np.loadtxt('tmp.txt', dtype=data, skiprows=True)
 
# Create a function
def gaussian(x, a, b, c):
    val = a * np.exp(-(x - b)**2 / c**2)
    return val
 
# Fit
popt, pcov = curve_fit(gaussian, enreg['distance'], enreg['effectif'])
 
# Print results
print "Scale =  %.3f +/- %.3f" % (popt[0], np.sqrt(pcov[0, 0]))
print "Offset = %.3f +/- %.3f" % (popt[1], np.sqrt(pcov[1, 1]))
print "Sigma =  %.3f +/- %.3f" % (popt[2], np.sqrt(pcov[2, 2]))
 
# Plot model
xm= np.linspace(enreg['distance'][0],enreg['distance'][-1], 100)
plt.plot(xm, gaussian(xm, popt[0], popt[1], popt[2]))
plt.show()

Je vais à présent essayer de le faire pour chaque pic, en modifiant ma variable xm, est ce comme cela que je dois procéder ?
En revanche, je ne pense pas avoir de bruid de fond dans mes données, du coup, je n'ai pas mis de sigma dans curve_fit ...

18/04/2013, 15h17
Isabella83

1 pièce(s) jointe(s)

En fait, j'avais des "trous" dans mes données, c'est à dire que quand je n'avais pas de distance, je ne mettais pas de 0.
Du coup, en complétant mon jeu de données, j'obtiens le graphiques mis en pièce jointe, qu'en pensez vous ?

1 pièce(s) jointe(s)

Salut,
En effet il n'y a pas de bruit de fond; 3 paramètres suffissent pour ajuster à une gaussienne.

J'ai corrigé ton programme, il manquait l'estimation des paramètres d'ajustement (argument p0 de curve_fit). Comme il y a plusieurs pics, la fonction ne sais pas lequel prendre.
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 #!/usr/bin/env python # -*- coding:Utf-8 -*- import numpy as np import matplotlib.pyplot as plt from scipy.optimize import curve_fit def gaussian(x, a, b, c): return a * np.exp(-(x - b)**2 / c**2) if __name__ == '__main__': data = np.dtype( [ ('distance',np.int), ('effectif',np.int) ] ) enreg = np.loadtxt('tmp.txt', dtype=data, skiprows=True) # Estimation des paramètres amplitude, position, largeur = 18, 4, 2 params = [ amplitude, position, largeur ] # Fit popt, pcov = curve_fit(gaussian, enreg['distance'], enreg['effectif'], p0=params) # Print results print "amplitude = %.3f +/- %.3f" % (popt[0], np.sqrt(pcov[0, 0])) print "position = %.3f +/- %.3f" % (popt[1], np.sqrt(pcov[1, 1])) print "largeur = %.3f +/- %.3f" % (popt[2], np.sqrt(pcov[2, 2])) # Plot model xm= np.linspace(enreg['distance'][0],enreg['distance'][-1], 512) plt.plot(xm, gaussian(xm, popt[0], popt[1], popt[2]), color='r', linewidth=2.0) plt.bar( enreg['distance'], enreg['effectif']) plt.show()
Il y a trop peu de statistique dans cette histogramme. Les pics ne ressemblent pas à des gaussiennes. Il faudrait plutôt essayer avec des classes de largeur 2. Ou alors les pics ont-ils vraiment cette forme ? Je n'arrive d'ailleurs pas à ajuster correctement les pics secondaires.

La méthode d'ajustement est meilleure que l'interpolation car elle élimine le bruit et les irrégularités non physiques, mais il y a trop peu de statistique pour que le résultat soit représentatif du signal d'origine. A toi de choisir.

19/04/2013, 10h20
Isabella83

Ah oui effectivement c'est mieux !
Est ce que je peux inclure plusieurs coordonnées de pic et ne faire qu'un seul schema ? Car j'ai une valeur de 0 après "l'amplitude" du pic, mais j'aimerai que la valeur de ne soit pas 0 ...