Führen Sie mit Python, matplotlib, einen logarithmischen Normalwahrscheinlichkeitsplot durch

Einführung

Der Bericht enthält die folgenden Daten.

Es scheint, dass die stochastische Flutungsrate durch eine logarithmische Normalverteilung mit zwei Variablen berechnet wird. Mit diesen Daten werden wir versuchen, Folgendes zu tun.

• 10000 Jahre Wahrscheinlichkeitsschätzung Hochwasserfluss
• Überprüfen Sie, ob es einer logarithmischen Normalverteilung mit zwei Variablen folgt

Es wird auf die Daten hochgerechnet, aber es ist besser, als die Zahlen in den begrenzten Informationen intuitiv anzugeben. Nun, es ist eine Programmpraxis, achten Sie also auf statistische Genauigkeit.

Ablauf der Prüfung

Vorbereitung für die Regressionsanalyse

Zur Vorbereitung der Regressionsanalyse berechnen wir die Prozentpunkte, die den nicht überschreitenden Wahrscheinlichkeiten in der Standardnormalverteilung entsprechen, und den gemeinsamen logarithmischen Wert der Durchflussrate.

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

Regressionsanalyse

Als nächstes wird eine lineare Regression mit der Zielvariablen als logarithmischem Wert der Durchflussrate und der erklärenden Variablen als Prozentpunkt der nicht übersteigenden Wahrscheinlichkeit durchgeführt, und die der 10.000-Jahres-Wahrscheinlichkeit entsprechende Hochwasserdurchflussrate wird geschätzt.

Berechnung des Regressionskoeffizienten und des Korrelationskoeffizienten durch lineare Regression mit 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)

Funktion für optimize.leastsq (ausgedrückt in einem Format gleich 0)

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

Schätzung der 1000-Jahres-Wahrscheinlichkeit anhand von Regressionsergebnissen

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

Zeichnen Sie auf Protokoll mit normaler Wahrscheinlichkeit

Ob es der logarithmischen Normalverteilung folgt oder nicht, wird auf einem logarithmischen Normalwahrscheinlichkeitspapier aufgezeichnet und bestätigt. Hier ist die horizontale Achse des Diagramms der Hochwasserfluss (logarithmischer Wert) und die vertikale Achse der Wahrscheinlichkeitswert. Daher wird die Standardachsenskala gelöscht und sowohl die horizontale als auch die vertikale Achse werden unabhängig voneinander festgelegt.

Löschen Sie Achsenskalen und Skalierungslinien

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

Das Diagramm wird als PNG-Bild gespeichert, die Ränder werden jedoch gelöscht, um es auf Qiita usw. einzufügen.

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

Programm

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()

Ausgabeergebnis

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

Referenzseite

Scipy statistische Funktionen http://kaisk.hatenadiary.com/entry/2015/02/17/192955

Linearer Regressionsfall von Scipy http://www2.kaiyodai.ac.jp/~kentaro/materials/new_HP/python/15fit_data3.html

Entfernen Sie die Bildränder mit matplotlib https://mzmttks.blogspot.my/2012/01/pylab-2.html