damyarou

python, GMT などのプログラム

設計 確率雨量の推定（１）計算式

設計

記事の最後に行く

はじめに

確率雨量の推定方法を紹介します。 この記述の元ネタは、以下の論文によります。

  • 星清：水文統計解析，開発土木研究所月報，No.540，1998年5月 （https://thesis.ceri.go.jp/db/files/0005005050.pdf）
  • 星清：現場のための水文統計(1)，開発土木研究所月報，No.540，1998年5月
  • 星清・新目竜一・宮原雅幸：現場のための水文統計(2)，開発土木研究所月報，No.541，1998年6月

ここでは、比較的良く使われると思われる、下記３種類の確率分布を考えます。 いづれも、３変数関数です。

  • ３変数対数正規分布 (Log-Normal distribution with 3 parameters, LN3)
  • 一般化極値分布（Generalized Extreme Value distribution, GEV)
  • 対数ピアソンIII型分布 (Log-Pearson type III distribution, LP3)

本文は、たまたま昔英文のものを作っていたので、それをコピーしています。 Jupyter の Markdown で正常に表示されることを確認して「はてな」に持ってきましたが、「はてな」では正常表示するのに苦労しました。 とりあえず、今後の参考のため、「はてな」の数式表示で気がついたことをメモしておきます。

  • 別行だての式では、{align}、{gather} が使える。 {equation} も使えるが {displaystyle} を使わないと行内表示となるためここでは使っていない。
  • 別行だて、行内表示とも、[ ..... ] は\[ ..... \]のように、必ずエスケープする。
  • 行内表示ではアンダーバー（_）は必要に応じてエスケープする。
  • 不等号は、うまく表示されない場合は、空白で囲む。
  • TeX記法での \text{ ..... }の中で [tex: ] は使えない。
  • 水圧鉄管の記事では、\text{ ... $...$ ... } は使えた。

以下、数式を羅列します。

Calculation Formulas

Basis

PDF and CDF

The relationship between PDF and CDF is shown below:

\begin{align} F(x)=\int_{-\infty}^x f(t) dt \end{align}


f(x) PDF (probability density function)
F(x) CDF (cumulative distribution function)

According to the definition of CDF, we can understand that F(x) : CDF is equal to non-exceedance probability for x.

Retuen period

For Annual Maximum Series data (AMS):

\begin{align} T=\frac{1}{1-p} \qquad p=1-\frac{1}{T} \end{align}

Plotting position formulas

A general formula for computing plotting position is defined as follow:

\begin{align} F[x_{(i)}] = \frac{ i - \alpha }{ N + 1 - 2 \alpha } \end{align}


N Number of sample
i Rank order in ascending order
x_{(i)} Value of sample with rank order i in ascending order
\alpha Constant for a paticular plotting position formula ( \alpha=0\sim 1)
F[x_{(i)}] Plotting position (like a non-exceedance probability)


Formula WeibullBlom CunnaneGringortenHazen
\alpha0 0.3750.40 0.44 0.5

Probability Distribution Model

Log-Normal Distribution with 3 parameters (LN3-distribution)

The parameter estimation method in the case of Hydrologic statistic $x$ following Log-Normal distribution with 3 parameters is shown below.

Probability Density Function

\begin{align} f(x)=\frac{1}{(x-a) \sigma_y \sqrt{2\pi}} \cdot \exp \left\{ -\frac{1}{2} \left[ \frac{\ln(x-a)-\mu_y}{\sigma_y}\right]^2\right\} \end{align}

Cumulative Distribution Function

\begin{align} F(x)=\Phi\left(\frac{\ln(x-a)-\mu_y}{\sigma_y}\right) \qquad \Phi(z)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^z\exp\left(-\frac{1}{2}t^2\right)dt \end{align}

Quantile x_p for non-exceedance probability p

\begin{align} z=\frac{\ln(x-a)-\mu_y}{\sigma_y} \qquad\rightarrow\qquad x=a+\exp(\mu_y+\sigma_y z) \end{align}
\begin{align} x_p=a+\exp(\mu_y+\sigma_y z_p) \qquad \text{$z_p$: value of $z$ at $p=\Phi(z)$} \end{align}

Estimation of parameters (Moment method)

\begin{align} \qquad \bar{x}=\frac{1}{N}\sum_{j=1}^N x_j \qquad S_x^2=\frac{1}{N}\sum_{j=1}^N (x_j-\bar{x})^2 \qquad C_{sx}=\frac{1}{N}\sum_{j=1}^N \left(\frac{x_j-\bar{x}}{S_x}\right)^3 \end{align}
\begin{align} \mu_x=\bar{x} \qquad \sigma_x=[N/(N-1)]^{1/2}S_x \qquad \gamma_x=\frac{\sqrt{N(N-1)}}{N-2}C_{sx} \end{align}

To get the unbiased skewness [tex\gamma_x], the following formula for Log-Normal distribution by Bobee and Robitaille is adopted. Note that C_{sx} to the power of 3 is multiplied by B.

\begin{gather} \gamma_x=C_{sx}(A+B\cdot C_{sx}{}^3) \\ \text{where,} \qquad A=1.01+7.01/N+14.66/N^2 \qquad B=1.69/N+74.66/N^2 \end{gather}

The relationship between moments and parameters is shown below:

\begin{align} &\mu_x=a+\exp(\mu_y)\cdot\exp(\sigma_y^2/2) \\ &\sigma_x=\exp(\mu_y)\sqrt{\exp(\sigma_y^2)\{\exp(\sigma_y^2)-1\}} \\ &\gamma_x=\{\exp(\sigma_y^2)+2\}\sqrt{\exp(\sigma_y^2)-1} \end{align}

The parameters which we want to know are \sigma_y, \mu_y and a in above equations. The third equation in above is expressed as follow,

\begin{align} \gamma_x=\{\exp(\sigma_y^2)+2\}\sqrt{\exp(\sigma_y^2)-1} \Rightarrow x^3+3x^2-4-\gamma_x^2=0 \quad \text{where,}\quad x=\exp(\sigma_y^2) \end{align}

The real positive solution for above cubic equation can be obtained by Cardano's method as follow,

\begin{align} x=\left(\beta+\sqrt{\beta^2-1}\right)^{1/3}+\left(\beta-\sqrt{\beta^2-1}\right)^{1/3}-1 \quad \text{where,}\quad \beta=1+\frac{\gamma_x^2}{2} \end{align}

As a result, parameters can be estimated using following equations with x which is a real positive silution of cubic equation shown above.

\begin{align} \begin{cases} \sigma_y=\sqrt{\ln(x)}\\ \mu_y=\ln\left(\cfrac{\sigma_x}{\sqrt{x(x-1)}}\right) \\ a=\mu_x-\exp(\mu_y)\cdot\exp(\sigma_y^2/2) \end{cases} \end{align}

Log-Pearson type III distribution (LP3-distribution)

The parameter estimation method in the case of Hydrologic statistic x following Log-Pearson type III distribution is shown below.

Probability Density Function

\begin{align} f(x)=\frac{1}{|a|\cdot\Gamma(b)\cdot x} \left(\frac{\ln x-c}{a}\right)^{b-1} \cdot \exp\left(-\frac{\ln x-c}{a}\right) \qquad a > 0 : \exp(c) < x < \infty \end{align}

\tex:c] is a location parameter, a is a scale parameter, b is a shape parameter.

Cumulative Distribution Function

\begin{align} F(x)=G\left(\frac{\ln x-c}{a}\right) \qquad G(w)=\frac{1}{\Gamma(b)}\int_0^w t^{b-1}\exp(-t)dt \qquad (a>0) \end{align}

Quantile x_p for non-exceedance probability p

\begin{align} w=\frac{\ln x-c}{a} \qquad\rightarrow\qquad x=\exp(c+a w) \end{align}
\begin{align} x_p=\exp(c+a w_p) \qquad \text{$w_p$: the value of $w$ at $p=G(w)$} \end{align}

Estimation of parameters (Moment method)

\begin{gather} y_j=\ln x_j \\ \bar{y}=\frac{1}{N}\sum_{j=1}^N y_j \qquad S_y^2=\frac{1}{N}\sum_{j=1}^N (y_j-\bar{y})^2 \qquad C_{sy}=\frac{1}{N}\sum_{j=1}^N \left(\frac{y_j-\bar{y}}{S_y}\right)^3 \end{gather}
\begin{align} \mu_y=\bar{y} \qquad \sigma_y=[N/(N-1)]^{1/2}S_y \qquad \gamma_y=\frac{\sqrt{N(N-1)}}{N-2}C_{sy} \end{align}

To get the unbiased skewness \gamma_x, the following formula for Pearson type III distribution by Bobee and Robitaille is adopted. Note that C_{sy} to the power of 2 is multiplied by B.

\begin{gather} \gamma_y=C_{sy}(A+B\cdot C_{sy}{}^2) \\ \text{where,} \qquad A=1+6.51/N+20.2/N^2 \qquad B=1.48/N+6.77/N^2 \end{gather}

The relationship between moments and parameters is shown below:

\begin{align} \mu_y=c+a \cdot b \qquad \sigma_y{}^2=a^2\cdot b \qquad \gamma_y=\frac{2 a}{|a|\sqrt{b}} \end{align}

As a result, parameters can be estimated using following equations.

\begin{align} \begin{cases} b=4/\gamma_y{}^2 \qquad (b>0) \\ a=\sigma_y/\sqrt{b} \qquad (\gamma_y<0 \rightarrow a=-\sigma_y/\sqrt{b}<0) \\ c=\mu_y-a b \end{cases} \end{align}

Note that if \gamma_y is less than 0, a shall be less than 0 and the value of w_p is for the value of 1-p.

Generalized Extreme Value distribution (GEV-distribution)

The parameter estimation method in the case of Hydrologic statistic x following Generalized Extreme Value distribution is shown below. Gumbel distribution is the same as Generalized Extreme value distribution with k=0.

Probability Density Function

\begin{align} f(x)=\frac{1}{a}\left(1-k\frac{x-c}{a}\right)^{1/k-1}\cdot\exp\left[-\left(1-k\frac{x-c}{a}\right)^{1/k}\right] \qquad (k\neq 0) \end{align}

c is a location parameter, a is a scale parameter, k is a shape parameter.

Cumulative Distribution Function

\begin{align} F(x)=\exp\left[-\left(1-k\frac{x-c}{a}\right)^{1/k}\right] \qquad (k\neq 0) \end{align}

Quantile x_p for non-exceedance probability p

\begin{align} p=\exp\left[-\left(1-k\frac{x-c}{a}\right)^{1/k}\right] \qquad\rightarrow\qquad x=c+\frac{a}{k}\cdot\left\{1-[-\ln(p)]^k\right\} \end{align}
\begin{align} x_p=c+\frac{a}{k}\cdot\left\{1-[-\ln(p)]^k\right\} \end{align}

Estimation of parameters (L-moments method)

\begin{align} & b_0=\frac{1}{N}\sum_{j=1}^N x_{(j)} \\ & b_1=\frac{1}{N(N-1)}\sum_{j=1}^N (j-1) x_{(j)} \\ & b_2=\frac{1}{N(N-1)(N-2)}\sum_{j=1}^N (j-1)(j-2) x_{(j)} \end{align}

where, x_{(i)} is a j-th value of sample x in ascending order. It means that x_{(1)} is the minimum value of sample data and x_{(N)} is the maximum value of sample data.

The values of \lambda_i can be calculated using following relationship by deeming b_i=\beta_i.

\begin{align} \lambda_1=\beta_0 \qquad \lambda_2=2\beta_1-\beta_0 \qquad \lambda_3=6\beta_2-6\beta_1+\beta_0 \end{align}

Parameters can be estimated using following relationship between L-Moments and parameters.

\begin{align} \begin{cases} k=7.8590d+2.9554d^2 \qquad \text{where, } \qquad d=\cfrac{2\lambda_2}{\lambda_3+3\lambda_2}-\cfrac{\ln(2)}{\ln(3)} \\ a=\cfrac{k\lambda_2}{(1-2^{-k})\cdot\Gamma(1+k)} \\ c=\lambda_1-\cfrac{a}{k}\cdot\left[1-\Gamma(1+k)\right] \end{cases} \end{align}

Thank you.

記事の先頭に行く