我已经计算了桥梁的荷载,我想用最大似然估计将Gumbel分布拟合到最高的20%。我需要帮助计算分布的参数。我已经读完了scipy.optimize公司但我无法理解如何应用其中的函数来估计双参数函数。在
以下是一些可能有帮助的理论: 有两个似然函数(L1和L2),一个用于高于某个阈值(x>;=C)的值,另一个用于低于(x<;C)的值,现在最有可能的参数是在两个函数max(L1*L2)之间乘法的最大值处的参数。在这种情况下,L1仍然是在席曦上概率密度函数值乘积的乘积,但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)))
首先,使用
scipy.optimize
优化函数的最简单方法是构造目标函数,以便第一个参数是需要优化的参数列表,而下面的参数指定其他内容,例如数据和固定参数。在其次,使用
numpy
提供的矢量化将非常有帮助因此,我们有这些:
在
trunc_GBL
函数中,我用缩放的pdf替换了您的pdf请看这里的基本原理,基本上是因为你的
L1
是基于pdf的,L2
是基于cdf的:http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_lifereg_sect018.htm然后我们注意到一个问题,请参阅最后一个输出中的
Current function value: 0.000000
。负对数似然函数为0。在这是因为:
^{2}$实际上是0。这意味着,根据您刚刚描述的模型,当阈值足够低时,总是会达到最大值,使得}为1(},对于数据中的所有项)。在
L1
不存在(x < threshold
为空)且{1-F(C)
为{因为这个原因,我觉得你的模特不太合适。你可能需要重新考虑一下。在
编辑
我们可以进一步分离
threshold
,并将其视为固定参数:并以不同的方式调用优化器:
这样,如果你想要70%的分位数,你可以简单地把它改成
np.percentile(X, 30)
等等。np.percentile()
只是另一种方法.quantile(0.8)
相关问题 更多 >
编程相关推荐