Python: Diagramme de l'état de distribution des données bidimensionnelles (donnant du sens aux contours de l'estimation de la densité du noyau)
introduction
Dans Article précédent, l'état de distribution des données bidimensionnelles a été tracé par contour basé sur l'estimation de la densité du noyau, mais la hauteur pour tracer la ligne de contour est particulièrement probable. Il a été décidé en fonction de la hauteur de la fonction de densité de probabilité sans tenir compte de la signification. En ce qui concerne le contour de la distribution de densité du noyau, je vois souvent des articles tels que "Si vous utilisez Seaborn, vous pouvez appliquer le contour de la distribution de données bidimensionnelle", mais j'ai vu celui qui expliquait ce que signifie la ligne. Jamais. De plus, si vous montrez les contours aux clients, vous pouvez voir qu'on leur demande: "Quelle est la signification technique de la ligne?" Par conséquent, cette fois, j'ai essayé de dessiner la ligne de contour en fonction de la valeur de probabilité en la développant un peu à partir de la fois précédente. Puisque je le fais avec ma propre compréhension, il peut y avoir des erreurs, mais dans ce cas, je voudrais signaler.
Comment faire
Vous trouverez ci-dessous le code de la partie où les contours sont soustraits de l'estimation de la densité du noyau. Dans l'estimation de la densité du noyau, une fonction de densité de probabilité est produite en fonction de la distribution des données (variable `` f '' sur la 6ème ligne du code ci-dessous).
01: xx,yy = np.mgrid[xmin:xmax:0.01,ymin:ymax:0.01]
02: positions = np.vstack([xx.ravel(),yy.ravel()])
03: value = np.vstack([x,y])
04: #kernel = gaussian_kde(value, bw_method='scott')
05: kernel = gaussian_kde(value, bw_method=0.1)
06: f = np.reshape(kernel(positions).T, xx.shape)
07: ll=fprop(f)
08: f=f/np.max(f)
09: plt.contourf(xx,yy,f,cmap=cm.Purples,levels=ll)
10: plt.contour(xx,yy,f,colors='#000080',linewidths=0.7,levels=ll)
Ici, la fonction fprop sur la 7ème ligne du code ci-dessus est introduite pour calculer la liste ll '' de la hauteur pour dessiner la ligne de contour. La fonction fpropest illustrée ci-dessous. Puisque nous utiliserons l'interpolation linéaire de scipy`` , importez scipy.interpolate```.
01: def fprop(f):
02: ff=np.ravel(f)
03: fmax=np.max(ff)
04: fsum=np.sum(ff)
05: asx=np.array([0.001,0.01,0.05,0.10,0.30,0.50,0.70,0.90,0.99])
06: asy=np.zeros(len(asx),dtype=np.float64)
07: for i,sx in enumerate(asx):
08: ii=np.where(sx<=ff/fmax); asy[i]=np.sum(ff[ii])/fsum
09: fip = interpolate.interp1d(asy, asx)
10: pr=np.array([0.99,0.90,0.70,0.50,0.30,0.10]) # probability
11: ch = fip(pr) # conter height corresponding to probability
12: con_h = ch.tolist()+[1.0] # np.array to list for plotting
13: print('fmax=',fmax) # peak value of probability density function
14: print('fsum=',fsum) # total sum of probability density function
15: print('pr=',pr) # specified probability
16: print('ch=',ch) # height of contour line corresponding to the specified probability
17: return con_h
La fonction ci-dessus fprop``` est une valeur de fonction de densité de probabilité` donnée aux coordonnées correspondant aux coordonnées bidimensionnelles spécifiées (xx, yy) (première ligne du code de dessin). A partir de f '', la hauteur (liste de valeurs de sortie: con_h```) à contourner en fonction de la valeur de probabilité donnée est spécifiée.
Ce que nous faisons est le suivant.
- 2ème ligne: Convertit le tableau 2D
- 3e ligne: Calculez la valeur maximale de la fonction de densité de probabilité (`` fmax```).
- 4ème ligne: Calculez la somme des valeurs de la fonction de densité de probabilité (
fsum```). Dans le cas de ce programme, la valeur sortie dans le programme estfsum = 9999.9999 ..., mais le maillage pour le calcul de la valeur de probabilité est donné par incréments de 0,01 (première ligne du code de dessin). ), La valeur intégrée de la fonction de densité de probabilité est donc9999,9999 x 0,01 x 0,01 = 0,9999 '', ce qui est presque 1, ce qui est logique. - 5ème ligne: spécifiez la hauteur pour dessiner le contour comme valeur temporaire.
- 6e à 8e lignes: calculez la probabilité que la zone soit inférieure à la hauteur de contour supposée.
- Lignes 9-11: Estimer la hauteur du contour (séquence
) correspondant à la probabilité non dépassée donnée (séquence` `` pr) par interpolation linéaire. Par conséquent, la relation entre la hauteur à laquelle les contours sont dessinés et la probabilité est, bien entendu, une valeur approximative. - 12ème ligne: ajoutez la hauteur maximale du contour pour le dessin.
En conséquence, le programme peut dessiner de manière probabiliste des lignes de contour avec des valeurs de 99%, 90%, 70%, 50%, 30% et 10% dans cette plage.
Exemple de dessin (bw_method = 0.1)
Lorsque le nombre de données est petit, la forme devient assez compliquée, mais le contour (ligne à 99%) correspond à peu près à la distribution des données et la figure est comme prévu. En passant, si cela fonctionne, cela ressemble à des nuages de prévisions météorologiques ou à des données Amedas.
programme
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from scipy.stats import gaussian_kde
from scipy import interpolate
def fprop(f):
ff=np.ravel(f)
fmax=np.max(ff)
fsum=np.sum(ff)
asx=np.array([0.001,0.01,0.05,0.10,0.30,0.50,0.70,0.90,0.99])
asy=np.zeros(len(asx),dtype=np.float64)
for i,sx in enumerate(asx):
ii=np.where(sx<=ff/fmax); asy[i]=np.sum(ff[ii])/fsum
fip = interpolate.interp1d(asy, asx)
pr=np.array([0.99,0.90,0.70,0.50,0.30,0.10]) # probability
ch = fip(pr) # conter height corresponding to probability
con_h = ch.tolist()+[1.0] # np.array to list for plotting
print('fmax=',fmax) # peak value of probability density function
print('fsum=',fsum) # total sum of probability density function
print('pr=',pr) # specified probability
print('ch=',ch) # height of contour line corresponding to the specified probability
return con_h
def sreq(x,y):
# y=a*x+b
res=np.polyfit(x, y, 1)
a=res[0]
b=res[1]
coef = np.corrcoef(x, y)
r=coef[0,1]
return a,b,r
def inp_obs():
fnameR='df_rgs1_tank_inp.csv'
df1_obs=pd.read_csv(fnameR, header=0, index_col=0) # read excel data
df1_obs.index = pd.to_datetime(df1_obs.index, format='%Y/%m/%d')
df1_obs=df1_obs['2016/01/01':'2018/12/31']
#
fnameR='df_rgs2_tank_inp.csv'
df2_obs=pd.read_csv(fnameR, header=0, index_col=0) # read excel data
df2_obs.index = pd.to_datetime(df2_obs.index, format='%Y/%m/%d')
df2_obs=df2_obs['2016/01/01':'2018/12/31']
#
return df1_obs,df2_obs
def inp_tank():
fnameR='df_rgs1_tank_result.csv'
df1_tank=pd.read_csv(fnameR, header=0, index_col=0) # read excel data
df1_tank.index = pd.to_datetime(df1_tank.index, format='%Y/%m/%d')
df1_tank=df1_tank['2001/01/01':'2018/12/31']
#
fnameR='df_rgs2_tank_result.csv'
df2_tank=pd.read_csv(fnameR, header=0, index_col=0) # read excel data
df2_tank.index = pd.to_datetime(df2_tank.index, format='%Y/%m/%d')
df2_tank=df2_tank['2001/01/01':'2018/12/31']
#
return df1_tank,df2_tank
def inp_mlp():
fnameR='df_rgs1_mlp_result.csv'
df1_mlp=pd.read_csv(fnameR, header=0, index_col=0) # read excel data
df1_mlp.index = pd.to_datetime(df1_mlp.index, format='%Y/%m/%d')
df1_mlp=df1_mlp['2001/01/01':'2018/12/31']
#
fnameR='df_rgs2_mlp_result.csv'
df2_mlp=pd.read_csv(fnameR, header=0, index_col=0) # read excel data
df2_mlp.index = pd.to_datetime(df2_mlp.index, format='%Y/%m/%d')
df2_mlp=df2_mlp['2001/01/01':'2018/12/31']
#
return df1_mlp,df2_mlp
def drawfig(x,y,sstr,fnameF):
fsz=12
xmin=0; xmax=3
ymin=0; ymax=3
plt.figure(figsize=(12,6),facecolor='w')
plt.rcParams['font.size']=fsz
plt.rcParams['font.family']='sans-serif'
#
for iii in [1,2]:
nplt=120+iii
plt.subplot(nplt)
plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
plt.xlabel('Discharge at RGS2 $Q_2$ (m$^3$/s)')
plt.ylabel('Discharge at RGS1 $Q_1$ (m$^3$/s)')
plt.gca().set_aspect('equal',adjustable='box')
plt.xticks([0,1,2,3], [1,10,100,1000])
plt.yticks([0,1,2,3], [1,10,100,1000])
bb=np.array([1,2,3,4,5,6,7,8,9,10,20,30,40,50,60,70,80,90,100,200,300,400,500,600,700,800,900,1000])
bl=np.log10(bb)
for a0 in bl:
plt.plot([xmin,xmax],[a0,a0],'-',color='#999999',lw=0.5)
plt.plot([a0,a0],[ymin,ymax],'-',color='#999999',lw=0.5)
plt.plot([xmin,xmax],[ymin,ymax],'-',color='#999999',lw=1)
#
if iii==1:
plt.plot(x,y,'o',ms=2,color='#000080',markerfacecolor='#ffffff',markeredgewidth=0.5)
a,b,r=sreq(x,y)
x1=np.min(x)-0.1; y1=a*x1+b
x2=np.max(x)+0.1; y2=a*x2+b
plt.plot([x1,x2],[y1,y2],'-',color='#ff0000',lw=2)
tstr=sstr+'\n$\log(Q_1)=a * \log(Q_2) + b$'
tstr=tstr+'\na={0:.3f}, b={1:.3f} (r={2:.3f})'.format(a,b,r)
#
if iii==2:
tstr=sstr+'\n(kernel density estimation)'
xx,yy = np.mgrid[xmin:xmax:0.01,ymin:ymax:0.01]
positions = np.vstack([xx.ravel(),yy.ravel()])
value = np.vstack([x,y])
#kernel = gaussian_kde(value, bw_method='scott')
kernel = gaussian_kde(value, bw_method=0.1)
f = np.reshape(kernel(positions).T, xx.shape)
ll=fprop(f)
f=f/np.max(f)
plt.contourf(xx,yy,f,cmap=cm.Purples,levels=ll)
plt.contour(xx,yy,f,colors='#000080',linewidths=0.7,levels=ll)
xs=xmin+0.03*(xmax-xmin)
ys=ymin+0.97*(ymax-ymin)
plt.text(xs, ys, tstr, ha='left', va='top', rotation=0, size=fsz,
bbox=dict(boxstyle='square,pad=0.2', fc='#ffffff', ec='#000000', lw=1))
#
plt.tight_layout()
plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
plt.show()
def main():
df1_obs, df2_obs =inp_obs()
df1_tank,df2_tank=inp_tank()
df1_mlp ,df2_mlp =inp_mlp()
fnameF='fig_cor_obs.png'
sstr='Observed discharge (2016-2018)'
xx=df2_obs['Q'].values
yy=df1_obs['Q'].values
xx=np.log10(xx)
yy=np.log10(yy)
drawfig(xx,yy,sstr,fnameF)
fnameF='fig_cor_tank.png'
sstr='Tank model (2001-2018)'
xx=df2_tank['Q_pred'].values
yy=df1_tank['Q_pred'].values
xx=np.log10(xx)
yy=np.log10(yy)
drawfig(xx,yy,sstr,fnameF)
fnameF='fig_cor_mlp.png'
sstr='MLP (2001-2018)'
xx=df2_mlp['Q_pred'].values
yy=df1_mlp['Q_pred'].values
xx=np.log10(xx)
yy=np.log10(yy)
drawfig(xx,yy,sstr,fnameF)
#==============
# Execution
#==============
if __name__ == '__main__': main()
c'est tout
Recommended Posts