<jupyter_start><jupyter_text>Data loading<jupyter_code>avg_income = pd.read_csv('../data/average_monthly_household_income.csv')
avg_income.shape
avg_income.head()<jupyter_output><empty_output><jupyter_text>L1 loss as goodness-of-fit measure<jupyter_code>def objective(c, scale):
    return -l1_loss(avg_income['2020'], c, scale)
    
    
def l1_loss(data, c, scale):
    data = np.asarray(data)
    
    # data must be a 1x10 vector, with the i-th element being 
    # the average income of the i-th decile of the population
    assert data.size == 10
    
    # scale down the data and parameter for robust numerical behaviour
    factor = 10 ** math.floor(math.log10(data.max()))
    
    decile_avgs = weibull_decile_averages(c, scale/factor) * factor

    # calculate and return L1 loss
    return np.linalg.norm(data-decile_avgs, 1)


def weibull_decile_averages(c, scale):
    # find the deciles of the Weibull
    deciles = weibull_min.ppf(np.linspace(0, 1, 11), c, scale=scale)
    
    # calculate and return the average for each decile using numerical integration
    return np.array([
        integrate.quad(
            lambda x: x * weibull_min.pdf(x, c, scale=scale), 
            deciles[i], 
            deciles[i+1]
        )[0] / .1
        for i in range(deciles.size-1)
    ])<jupyter_output><empty_output><jupyter_text>Parameter optimization<jupyter_code>pbounds = {'c': (1.3, 1.6), 'scale': (8000, 15000)}

optimizer = BayesianOptimization(
    f=objective,
    pbounds=pbounds,
    random_state=1,
)

optimizer.maximize(init_points=50, n_iter=100)
optimizer.max
plt.bar(range(10), avg_income['2020'], alpha=0.5)
plt.bar(range(10), 
        weibull_decile_averages(optimizer.max['params']['c'], scale=optimizer.max['params']['scale']),
        alpha=0.5)
plt.legend(['data', 'model'])
plt.grid(axis='y')
plt.show()<jupyter_output><empty_output><jupyter_text>Income at 95% percentile<jupyter_code>weibull_min.ppf(.95, optimizer.max['params']['c'], scale=optimizer.max['params']['scale'])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13683178028445317, 'fraction_numerical': 0.029916625796959292, 'mean_word_length': 4.112781954887218, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '146868', 'license': 'MIT', 'n_tokens_mistral': 788, 'n_tokens_neox': 732, 'n_words': 145, 'path': 'notebooks/weibull_model.ipynb', 'repo_name': 'foxwood/sg-income-dist', 'size': '2039'}