2

橋の負荷を計算しました。最尤推定を使用して、ガンベルの分布をそれらの最大 20% に適合させたいと考えています。分布のパラメータを計算するのに助けが必要です。scipy.optimize のドキュメントを読みましたが、2 つのパラメーター関数を推定するためにそこに関数を適用する方法を理解できません。

役立つかもしれない理論は次のとおりです。2 つの尤度関数 (L1 と L2) があり、1 つはあるしきい値を超える値 (x>=C) 用で、もう 1 つはそれ以下の値 (x < C) 用です。 2 つの関数 max(L1*L2) 間の乗算の最大値。この場合、L1 はやはり xi での確率密度関数の値の乗算の積ですが、L2 はしきい値 C を超える確率 (1-F(C)) です。

ここに私が書いたいくつかのコードがあります:

non_truncated_data = ([15.999737471905252, 16.105716234887431, 17.947809230275304, 16.147752064149291, 15.991427126788327, 16.687542227378565, 17.125139229445359, 19.39645340792385, 16.837044960487795, 15.804473320190725, 16.018569387471025, 16.600876724289019, 16.161306985203151, 17.338636901595873, 18.477371969176406, 17.897236722220281, 16.626465201654593, 16.196548622931672, 16.013794215070927, 16.30367884232831, 17.182106070966608, 18.984566931768452, 16.885737663740024, 16.088051117522948, 15.790480003140173, 18.160947973898388, 18.318158853376037])

threshold = 15.78581825859324

def maximum_likelihood_function(non_truncated_loads, threshold, loc, scale):
    """Calculates maximum likelihood function's value for given truncated data
    with given parameters.
    Maximum likelihood function for truncated data is L1 * L2. Where L1 is a 
    product of multiplication of pdf values at non-truncated known values
    (non_truncated_values). L2 is a the probability that threshold value will
    be exceeded.
    """
        is_first = True
    # calculates L1
    for x in non_truncated_loads:
        if is_first:
            L1 = gumbel_pdf(x, loc, scale)
            is_first = False
        else:
            L1 *= gumbel_pdf(x, loc, scale)

    # calculates L2
    cdf_at_threshold = gumbel_cdf(threshold, loc, scale)
    L2 = 1 - cdf_at_threshold

    return L1*L2

def gumbel_pdf(x, loc, scale):
    """Returns the value of Gumbel's pdf with parameters loc and scale at x .
    """
    # exponent
    e = math.exp(1)
    # substitute
    z = (x - loc)/scale

    return (1/scale) * (e**(-(z + (e**(-z)))))

def gumbel_cdf(x, loc, scale):
    """Returns the value of Gumbel's cdf with parameters loc and scale at x.
    """
    # exponent
    e = math.exp(1)

    return (e**(-e**(-(x-loc)/scale)))
4

1 に答える 1

1

まず、 を使用して関数を最適化する最も簡単な方法はscipy.optimize、最初の引数が最適化する必要があるパラメーターのリストであり、次の引数がデータや固定パラメーターなどの他のものを指定するようにターゲット関数を作成することです。

第二に、によって提供されるベクトル化を使用すると非常に役立ちます。numpy

したがって、次のようになります。

In [61]:
#modified pdf and cdf
def gumbel_pdf(x, loc, scale):
    """Returns the value of Gumbel's pdf with parameters loc and scale at x .
    """
    # substitute
    z = (x - loc)/scale
 
    return (1./scale) * (np.exp(-(z + (np.exp(-z)))))
 
def gumbel_cdf(x, loc, scale):
    """Returns the value of Gumbel's cdf with parameters loc and scale at x.
    """
    return np.exp(-np.exp(-(x-loc)/scale))
In [62]:

def trunc_GBL(p, x):
    threshold=p[0]
    loc=p[1]
    scale=p[2]
    x1=x[x<threshold]
    nx2=len(x[x>=threshold])
    L1=(-np.log((gumbel_pdf(x1, loc, scale)/scale))).sum()
    L2=(-np.log(1-gumbel_cdf(threshold, loc, scale)))*nx2
    #print x1, nx2, L1, L2
    return L1+L2
In [63]:

import scipy.optimize as so
In [64]:
#first we make a simple Gumbel fit
so.fmin(lambda p, x: (-np.log(gumbel_pdf(x, p[0], p[1]))).sum(), [0.5,0.5], args=(np.array(non_truncated_data),))
Optimization terminated successfully.
         Current function value: 35.401255
         Iterations: 70
         Function evaluations: 133
Out[64]:
array([ 16.47028986,   0.72449091])
In [65]:
#then we use the result as starting value for your truncated Gumbel fit
so.fmin(trunc_GBL, [17, 16.47028986,   0.72449091],  args=(np.array(non_truncated_data),))
Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 25
         Function evaluations: 94
Out[65]:
array([ 13.41111111,  16.65329308,   0.79694   ])

ここでtrunc_GBL関数で、pdfをスケーリングされたpdfに置き換えました

ここに画像の説明を入力

ここで理論的根拠を参照してください。基本的には、あなたL1がpdfベースでL2cdfベースであるためです。

Current function value: 0.000000次に、最後の出力に問題があることに気付きます。負の対数尤度関数は 0 です。

それの訳は：

In [66]:

gumbel_cdf(13.41111111,  16.47028986,   0.72449091)
Out[66]:
2.3923515777163676e-30

L1事実上 0 です。これは、先ほど説明したモデルに基づいて、存在しない (x < thresholdが空である) およびL21 (データ内のすべてのアイテム1-F(C)がである)など、しきい値が十分に低い場合に常に最大値に達することを意味します。1

このため、あなたのモデルは私にはあまり正しく見えません。考え直した方がいいかもしれません。

編集

さらに分離thresholdして、固定パラメーターとして扱うことができます。

def trunc_GBL(p, x, threshold):
    loc=p[0]
    scale=p[1]
    x1=x[x<threshold]
    nx2=len(x[x>=threshold])
    L1=(-np.log((gumbel_pdf(x1, loc, scale)/scale))).sum()
    L2=(-np.log(1-gumbel_cdf(threshold, loc, scale)))*nx2
    #print x1, nx2, L1, L2
    return L1+L2

オプティマイザーを別の方法で呼び出します。

so.fmin(trunc_GBL, [0.5, 0.5], args=(X, np.percentile(X, 20)))
Optimization terminated successfully.
         Current function value: 20.412818
         Iterations: 72
         Function evaluations: 136
Out[9]:
array([ 16.34594943,   0.45253201])

このようにして、70% の分位点が必要な場合は、単純に変更することができますnp.percentile(X, 30)。np.percentile()別の方法です.quantile(0.8)

于 2014-04-22T21:06:52.787 に答える