免赔额建模困难使用scipy.optimize.minimize实现收敛

Question

我是Python新手。如果我的代码有问题，我会提前道歉。在过去的几周里，我一直在研究一个叫做免赔额建模的难题。我将尽量避免不必要的统计术语，并保持问题的简单，因为据我所知，我的问题是编程/优化问题，而不是统计问题

如果您认为应该，请将此线程移动到更合适的堆栈交换站点

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

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

我已经尝试过使用Nelder-Mead和SLSQP方法（注释掉了SLSQP），但是返回的值是笨拙和无意义的

以下是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估计量）获得的。我觉得我肯定还有别的地方做错了。而且，尽管听起来很蹩脚，我甚至使用Excel的解算器重现了完全相同的问题，并遇到了类似的数值收敛问题。我很乐意提供有关这个问题的补充细节，无论是技术性的还是非技术性的

Answer 1

Scipy包含许多优化算法，如果没有更好的建议，为什么不试试呢？你不能用解析的方法计算梯度这一事实是无关紧要的，scipy会用数值的方法来计算

具体来说，您可能希望尝试使用BFGS，如果这也失败了，那么您的代码中可能存在错误（我们很难发现，因为它目前无法运行），或者您的目标函数更难实现，需要一个更强大（尽管速度较慢）的全局优化器。Scipy附带了很多（差异进化、盆地跳跃、SHGO、双重退火），web上充满了Python的其他可选全局优化例程