AirTerminalDevice/model/utils/fit_utils.py

import numpy as np
import pandas as pd
import pymc as pm
from sklearn.metrics import r2_score,mean_absolute_error,mean_absolute_percentage_error

def record(name,var):
    pm.Deterministic(f'{name}_mu',var)

def observe(name,var,observed,sigma=1):
    if isinstance(observed,pd.DataFrame):
        observed = observed.loc[:,name].values
    mu    = pm.Deterministic(f'{name}_mu',var)
    sigma = pm.HalfNormal(f'{name}_sigma',sigma=sigma)
    pm.Normal(name,mu=mu,sigma=sigma,observed=observed)

def reorder_posterior(prior:dict,posterior:dict):
    param_posterior_reorder = {'F_air':{}}
    for equp_name in prior.keys():
        param_posterior_reorder.setdefault(equp_name,{})
        for param_name,param_value in posterior.items():
            if '__' in param_name:
                continue
            if param_name == 'F_air_val_rw':
                param_value = np.median(param_value[-5:])
            if param_name.startswith(equp_name):
                param_name_adj = param_name.replace(f'{equp_name}_','')
                param_posterior_reorder[equp_name][param_name_adj] = param_value
    return param_posterior_reorder

def get_fitted_result(
    posterior    : dict,
    observed_data: pd.DataFrame,
    plot_TVP     : bool
) -> tuple:
    # 样本内预测数据
    TVP_data = []
    for param_name in posterior.keys():
        if param_name.replace('_mu','') not in observed_data.columns:
            continue
        TVP_data.append(
            pd.DataFrame(
                {
                    'param_name': param_name.replace('_mu',''),
                    'real'      : observed_data.loc[:,param_name.replace('_mu','')].values,
                    'pred'      : posterior[param_name]
                }
            )
        )
    TVP_data = pd.concat(TVP_data,axis=0)
    
    group_by_data = TVP_data.groupby(['param_name'])[['pred','real']]
    TVP_metric = (
        pd.concat(
            [
                group_by_data.apply(lambda dt:r2_score(dt.real,dt.pred)),
                group_by_data.apply(lambda dt:mean_absolute_error(dt.real,dt.pred)),
                group_by_data.apply(lambda dt:mean_absolute_percentage_error(dt.real,dt.pred)),
            ],
            axis=1
        )
        .set_axis(['R2','MAE','MAPE'],axis=1)
        .sort_values(by='R2',ascending=True)
    )
    
    if plot_TVP:
        import plotnine as gg
        plot = (
            TVP_data
            .pipe(gg.ggplot)
            + gg.aes(x='real',y='pred')
            + gg.geom_point()
            + gg.facet_wrap(facets='param_name',scales='free')
            + gg.geom_abline(intercept=0,slope=1,color='red')
            + gg.theme(figure_size=[10,10])
        )
        plot.show()
        
    return TVP_data,TVP_metric