0

我是 Python 新手。如果我的代码感觉不舒服，我提前道歉。在过去的几周里，我一直在研究一项称为可扣除建模的困难统计挑战。我会尽量避免不必要的统计术语并保持问题简单，因为据我了解，我的问题是编程/优化问题，而不是统计问题。

如果您认为应该，请将此主题移至更合适的堆栈交换站点

我有 8 个参数，其中 2 个是受约束的，应该是正数（phi_freq和phi_sev）。我正在尝试最大化对数似然函数，该函数基本上是所有这些参数的非凸、多模态、连续、实值、不可微 (AFAIK) 函数。哎呀！

我知道这些问题意味着提供给搜索算法的种子值将对我们收敛到的局部/全局最优值产生巨大影响，但据我所知，我的起始值是可靠的，即使是硬编码的（在提供的示例中下面），是辅助优化的结果，而不仅仅是主观预感。

我曾尝试使用库的Nelder-MeadandSLSQP方法（SLSQP已注释掉）scipy.optimize.minimize，但返回的值很尴尬且毫无意义。

以下是 MWE：

from scipy.stats import gamma
from scipy import special
import pandas as pd
import numpy as np
import scipy as s

def freq_mean(intercept_freq, beta_veh_age_log, beta_dr_section_b):
    return np.exp(intercept_freq + beta_veh_age_log * freq_data['veh_age_log'] + beta_dr_section_b * freq_data['dr_section_b'])

def sev_mean(intercept_sev, alpha_veh_age_log, alpha_acv_log, df):
    return np.exp(intercept_sev + alpha_veh_age_log * df['veh_age_log'] + alpha_acv_log * df['acv_log'])

def freq_dist(x, mu, phi):
    return (sp.special.gamma(phi + x) / sp.special.gamma(phi) / sp.special.factorial(x)) * ((phi  / (phi + mu)) ** phi) * ((mu / (phi + mu)) ** x)

def sev_dist(x, mu, phi, ded):
    gamma_pdf = (((phi / mu) ** phi) / sp.special.gamma(phi)) * (x ** (phi - 1.0)) * np.exp(-phi * x / mu)
    gamma_sdf = 1.0 - sp.stats.gamma.cdf(x = phi * ded / mu, a = phi, scale = 1.0)

    if gamma_sdf == 0.0:
        return 0.0
    else:
        return gamma_pdf / gamma_sdf

def sev_loglik(intercept_sev, alpha_veh_age_log, alpha_acv_log, phi_sev):
    sev_data['mu_sev'] = sev_mean(intercept_sev, alpha_veh_age_log, alpha_acv_log, sev_data)
    return sev_data.apply(lambda x : np.log(sev_dist(x['claim_amount'], x['mu_sev'], phi_sev, x['deductible'])), axis = 1).sum()

def freq_loglik(intercept_freq, beta_veh_age_log, beta_dr_section_b, phi_freq, intercept_sev, alpha_veh_age_log, alpha_acv_log, phi_sev):
    freq_data['mu_sev'] = sev_mean(intercept_sev, alpha_veh_age_log, alpha_acv_log, freq_data)
    freq_data['claim_prob'] = 1.0 - sp.stats.gamma.cdf(x = phi_sev * freq_data['deductible'] / freq_data['mu_sev'], a = phi_sev, scale = 1.0)
    freq_data['mu_freq'] = freq_mean(intercept_sev, beta_veh_age_log, beta_dr_section_b)
    return freq_data.apply(lambda x : np.log(freq_dist(x['claim_count'], x['exposure'] * x['claim_prob'] * x['mu_freq'], x['exposure'] * phi_freq)), axis = 1).sum()

def obj_func(param):
    intercept_freq, beta_veh_age_log, beta_dr_section_b, phi_freq, intercept_sev, alpha_veh_age_log, alpha_acv_log, phi_sev = param
    ll_sev = sev_loglik(intercept_sev, alpha_veh_age_log, alpha_acv_log, phi_sev)
    ll_freq = freq_loglik(intercept_freq, beta_veh_age_log, beta_dr_section_b, phi_freq, intercept_sev, alpha_veh_age_log, alpha_acv_log, phi_sev)
    return -(ll_freq + ll_sev)

def mle(nfev = 100):    
    intercept_freq, beta_veh_age_log, beta_dr_section_b, phi_freq, intercept_sev, alpha_veh_age_log, alpha_acv_log, phi_sev = (-1.87515443, -0.5675389200, -0.0216802900, 23.78568667, 10.42040743, 0.00465891, 0.00072216, 0.69497075)
    seeds = np.array([intercept_freq, beta_veh_age_log, beta_dr_section_b, phi_freq, intercept_sev, alpha_veh_age_log, alpha_acv_log, phi_sev])
    return sp.optimize.minimize(fun = obj_func, x0 = seeds, method = 'Nelder-Mead', options = {'maxfev' : nfev})
    #cons = ({'type': 'ineq', 'fun' : lambda x : x[3]}, {'type': 'ineq', 'fun' : lambda x : x[7]})
    #return sp.optimize.minimize(fun = obj_func, x0 = seeds, method = 'SLSQP', constraints = cons)

policies = pd.read_csv('C:/policies.txt', sep = '\t')
claims = pd.read_csv('C:/claims.txt', sep = '\t')

freq_data = pd.merge(left = policies, right = claims.groupby('ID').agg(claim_count = ('claim_amount', 'count')), how = 'left', on = 'ID')
freq_data['claim_count'].fillna(0, inplace = True)
sev_data = pd.merge(left = claims[['ID', 'claim_amount']], right = policies.drop(['dr_section_b', 'exposure'], axis = 1), how = 'left', on = 'ID')

opt = mle(4000)

我不确定目前如何以及是否需要提供数据（示例代码中引用的policies和claims平面文件）。我准备这样做，但需要先将它们匿名化。

所以，我想此时任何指针都会受到欢迎，因为我快要碰壁了。有一个全局解决方案，必须有，这意味着mle（最大似然估计）必须存在。我的种子值非常现实，因为它们是通过矩匹配（所谓的method of moments估计器）获得的。我觉得一定有什么我做错了。此外，尽管听起来很蹩脚，但我什至使用 's 求解器重现了完全相同的问题Excel，并遇到了类似的数值收敛问题。我很乐意提供有关此问题的补充细节，无论技术与否。

4

1 回答 1

1

Scipy 包含许多优化算法，由于缺乏更好的建议，您为什么不尝试它们呢？您无法解析地计算梯度这一事实是无关紧要的，scipy 会在数值上为您做到这一点。

具体来说，您可能想尝试使用 BFGS，如果这也失败了，那么您的代码中可能存在错误（我们无法轻易找到它，因为它目前无法运行）或者您的目标函数更难它需要一个更强大——虽然速度更慢——的全局优化器。Scipy 附带了很多（差分进化、盆地跳跃、SHGO、双重退火），并且网络上充满了 Python 的其他替代全局优化例程。

于 2019-12-18T06:30:13.573 回答