Effectuer un tracé de probabilité normale logarithmique avec Python, matplotlib

introduction

Le rapport contient les données suivantes.

Il semble que le débit de crue stochastique soit calculé par une distribution normale logarithmique à deux variables. En utilisant ces données, nous essaierons de faire ce qui suit.

• 10000 ans d'écoulement de crue
• Vérifiez s'il suit une distribution normale logarithmique à deux variables

Il sera extrapolé aux données, mais c'est mieux que de donner les chiffres par intuition dans les informations limitées. Eh bien, c'est une pratique de programme, alors gardez un œil sur la rigueur statistique.

Déroulement de l'examen

Préparation à l'analyse de régression

En préparation de l'analyse de régression, les points de pourcentage correspondant aux probabilités non dépassantes dans la distribution normale standard et la valeur logarithmique commune du débit sont calculés.

ppp=norm.ppf(Pin, loc=0, scale=1)
qqq=np.log10(Qin)

analyse de régression

Ensuite, une régression linéaire est effectuée avec la variable objective comme valeur logarithmique du débit et la variable explicative comme point de pourcentage de la probabilité de non-dépassement, et le débit de crue correspondant à la probabilité de 10000 ans est estimé.

Calcul du coefficient de régression et du coefficient de corrélation par régression linéaire à l'aide de scipy: Optimize.leastsq

parameter0 = [0.0,0.0]
result = optimize.leastsq(func,parameter0,args=(ppp,qqq))
aa=result[0][0]
bb=result[0][1]
rr=np.corrcoef(ppp,qqq)

Fonction pour Optimize.leastsq (exprimé dans un format égal à 0)

def func(parameter,x,y):
    a = parameter[0]
    b = parameter[1]
    residual = y-(a*x+b)
    return residual

Estimation de la probabilité sur 1000 ans à l'aide des résultats de la régression

pp=0.9999
xq=norm.ppf(pp, loc=0, scale=1)
Qpp=10**(aa*xq+bb)

Tracé sur papier de probabilité log normale

Si elle suit ou non la distribution log normale est tracée sur un papier de probabilité log normale et confirmée. Ici, l'axe horizontal du graphique est le débit de crue (valeur logarithmique) et l'axe vertical est la valeur de probabilité. Par conséquent, l'échelle de l'axe par défaut est supprimée et l'axe horizontal et l'axe vertical sont définis indépendamment.

Supprimer les échelles des axes et les lignes d'échelle

plt.tick_params(labelbottom='off')
plt.tick_params(labelleft='off')
plt.tick_params(which='both', width=0)

Le graphique est enregistré en tant qu'image png, mais les marges sont supprimées en cas de collage sur Qiita, etc.

plt.savefig('fig_flood_p.png', bbox_inches="tight", pad_inches=0.2)

programme

py_qq.py

import numpy as np
from scipy.stats import norm # normal distribution
from scipy import optimize
import matplotlib.pyplot as plt


# function for optimize.leastsq
def func(parameter,x,y):
    a = parameter[0]
    b = parameter[1]
    residual = y-(a*x+b)
    return residual

#==============================
# data analysis
#==============================
# flood discharge data
Qin=np.array([950,1650,2190,2780,3610,4330,5090,6190,7090])
# non-exceedance probability data
Pin=np.array([0.50,0.80,0.90,0.95,0.98,0.99,0.995,0.998,0.999])

ppp=norm.ppf(Pin, loc=0, scale=1)
qqq=np.log10(Qin)

# least square method
parameter0 = [0.0,0.0]
result = optimize.leastsq(func,parameter0,args=(ppp,qqq))
aa=result[0][0]
bb=result[0][1]
rr=np.corrcoef(ppp,qqq)

# estimation of 10,000 years return period flood
pp=0.9999
xq=norm.ppf(pp, loc=0, scale=1)
Qpp=10**(aa*xq+bb)

print('log10(Q) = a * x + b')
print('a={aa:10.6f} b={bb:10.6f} r={rr[0][1]:10.6f}'.format(**locals()))
print('Q={Qpp:10.3f} (pp={pp:0.4f})'.format(**locals()))


#==============================
# plot by matplotlib
#==============================
fig = plt.figure(figsize=(7,8))
xmin=np.log10(100)
xmax=np.log10(20000)
ymin=norm.ppf(0.0001, loc=0, scale=1)
ymax=norm.ppf(0.9999, loc=0, scale=1)
plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
plt.tick_params(labelbottom='off')
plt.tick_params(labelleft='off')
plt.tick_params(which='both', width=0)

# data plot
plt.plot(qqq,ppp,'o')
# straight line by regression analysis
plt.plot([xmin, xmax], [(xmin-bb)/aa, (xmax-bb)/aa], color='k', linestyle='-', linewidth=1)

# setting of x-y axes
_dy=np.array([0.0001,0.001,0.01,0.1,0.5,0.9,0.99,0.999,0.9999])
dy=norm.ppf(_dy, loc=0, scale=1)
plt.hlines(dy, xmin, xmax, color='grey')
_dx=np.array([100,200,300,400,500,600,700,800,900,1000,2000,3000,4000,5000,6000,7000,8000,9000,10000,20000])
dx=np.log10(_dx)
plt.vlines(dx, ymin, ymax, color='grey')
fs=16
for i in range(2,5):
    plt.text(float(i), ymin-0.1, str(10**i), ha = 'center', va = 'top', fontsize=fs)
for i in range(0,9):
    plt.text(xmin-0.01, dy[i], str(_dy[i]), ha = 'right', va = 'center', fontsize=fs)
plt.text(0.5*(xmin+xmax), ymin-0.5, 'Flood discharge (m$^3$/s)', ha = 'center', va = 'center', fontsize=fs)
plt.text(xmin-0.25,0.5*(ymin+ymax),'Non-exceedance probability', ha = 'center', va = 'center', fontsize=fs, rotation=90)

# image saving and image showing
plt.savefig('fig_flood_p.png', bbox_inches="tight", pad_inches=0.2)
plt.show()

Résultat de sortie

$ python3 py_qq.py
log10(Q) = a * x + b
a=  0.282460 b=  2.978647 r=  0.999996
Q= 10693.521 (pp=0.9999)

Site de référence

Fonctions statistiques Scipy http://kaisk.hatenadiary.com/entry/2015/02/17/192955

Cas de régression linéaire par Scipy http://www2.kaiyodai.ac.jp/~kentaro/materials/new_HP/python/15fit_data3.html

Supprimer les marges de l'image avec matplotlib https://mzmttks.blogspot.my/2012/01/pylab-2.html