0

私はPythonが初めてです。私のコードがぎこちなく感じる場合は、事前にお詫び申し上げます。ここ数週間、私は控除可能なモデリングと呼ばれる難しい統計上の課題に取り組んできました。私が理解しているように、私の問題はプログラミング/最適化の問題であり、統計の問題ではありません。

必要に応じて、このスレッドをより適切なスタック交換サイトに移動してください。

8 つのパラメーターがあり、そのうちの 2 つは制約されており、正 (phi_freqおよびphi_sev) である必要があります。基本的に、これらすべてのパラメーターの非凸、マルチモーダル、連続、実数値、非微分 (AFAIK) 関数である対数尤度関数を最大化しようとしています。うわぁ！

このような問題は、検索アルゴリズムに提供されるシード値が、収束するローカル/グローバル最適値に多大な影響を与えることを意味しますが、私の知る限り、私の開始値はしっかりしており、ハードコーディングされています (提供されたサンプルでは)以下)、単なる主観的な直感ではなく、補助的な最適化の結果でした。

Nelder-MeadライブラリのandSLSQPメソッド (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推定値) によって取得されたため、非常に現実的です。私が間違っていることが他にあるに違いないと感じています。また、不自由に聞こえるかもしれませんが、Excelのソルバーを使用してまったく同じ問題を再現し、同様の数値収束の問題に遭遇しました。この問題に関して、技術的かどうかにかかわらず、補足的な詳細を提供させていただきます。

4

1 に答える 1

1

Scipy には多数の最適化アルゴリズムが含まれています。より良い提案がないため、それらを試してみませんか? 勾配を分析的に計算できないという事実は関係ありません.scipyは数値的にそれを行います.

具体的には、BFGS を試して使用することをお勧めします。これも失敗する場合は、コードに誤りがあるか (このままでは実行できないため、簡単には見つけることができません)、目的関数がはるかに難しいかのいずれかです。また、より強力な (低速ではありますが) グローバル オプティマイザが必要です。Scipy にはそれらのかなりの数 (差分進化、ベイスン ホッピング、SHGO、デュアル アニーリング) が付属しており、Web には Python 用の他の代替グローバル最適化ルーチンがたくさんあります。

于 2019-12-18T06:30:13.573 に答える