Aula 13 - Análise de Sensibilidade.ipynb

Kernel: SageMath 9.5

Análise de Sensibilidade

Este tipo de análise se propõe a medir a influencia de cada parâmetro do modelo sobre o resultado da simulação. A Sensibilidade do modelo à cada parâmetros é representada por um valor numérico chamado de índice de sensibilidade. estes índices pode ser de vários tipos:

• Indices de primeira ordem medem a contribuição de cada parâmetro individualmente para a variância da saída do modelo.$$S_i=\frac{V_i}{Var(Y)}$$

• Indices de segunda ordem medem a contribuição de pares de parâmetros para a variância da saída do modelo.

• Indices de ordem total medem a contribuição de cada parâmetro incluindo todas as suas interações para a variância da saída do modelo. $${\displaystyle S_{Ti}={\frac {E_{{\textbf {X}}_{\sim i}}\left(\operatorname {Var} _{X_{i}}(Y\mid \mathbf {X} _{\sim i})\right)}{\operatorname {Var} (Y)}}=1-{\frac {\operatorname {Var} _{{\textbf {X}}_{\sim i}}\left(E_{X_{i}}(Y\mid \mathbf {X} _{\sim i})\right)}{\operatorname {Var} (Y)}}}$$

Para realizarmos esta análise vamos precisar instalar duas bibliotecas: salib e epimodels. Para instalá-los, basta ativar a shell do Sage:

sage -sh

e executar os seguintes comandos:

pip install salib
pip install epimodels
pip install seaborn
exit

Ou diretamente aqui no notebook:

!pip3 install -U salib epimodels==0.4.0 seaborn

Defaulting to user installation because normal site-packages is not writeable
Requirement already satisfied: salib in /home/fccoelho/.sage/local/lib/python3.10/site-packages (1.4.6)
Collecting epimodels==0.4.0
  Downloading epimodels-0.4.0-py3-none-any.whl (10 kB)
Requirement already satisfied: seaborn in /home/fccoelho/.sage/local/lib/python3.10/site-packages (0.12.0)
Collecting matplotlib<4.0.0,>=3.6.1
  Using cached matplotlib-3.6.1-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (11.8 MB)
Collecting mypy==0.982
  Using cached mypy-0.982-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (17.7 MB)
Requirement already satisfied: sphinx in /usr/lib/python3/dist-packages (from epimodels==0.4.0) (4.3.2)
Collecting sympy<2.0.0,>=1.11.1
  Using cached sympy-1.11.1-py3-none-any.whl (6.5 MB)
Collecting cython==0.29.32
  Using cached Cython-0.29.32-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.manylinux_2_24_x86_64.whl (1.9 MB)
Collecting numpy<2.0.0,>=1.23.3
  Using cached numpy-1.23.3-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (17.1 MB)
Collecting scipy<2.0.0,>=1.9.2
  Using cached scipy-1.9.2-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (33.7 MB)
Collecting pyitlib
  Using cached pyitlib-0.2.2-py3-none-any.whl
Requirement already satisfied: typing-extensions>=3.10 in /usr/local/lib/python3.10/dist-packages (from mypy==0.982->epimodels==0.4.0) (4.3.0)
Requirement already satisfied: tomli>=1.1.0 in /usr/local/lib/python3.10/dist-packages (from mypy==0.982->epimodels==0.4.0) (2.0.1)
Requirement already satisfied: mypy-extensions>=0.4.3 in /usr/local/lib/python3.10/dist-packages (from mypy==0.982->epimodels==0.4.0) (0.4.3)
Requirement already satisfied: multiprocess in /home/fccoelho/.sage/local/lib/python3.10/site-packages (from salib) (0.70.13)
Requirement already satisfied: pandas>=1.1.2 in /usr/local/lib/python3.10/dist-packages (from salib) (1.4.3)
Requirement already satisfied: packaging>=20.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib<4.0.0,>=3.6.1->epimodels==0.4.0) (21.3)
Requirement already satisfied: cycler>=0.10 in /usr/lib/python3/dist-packages (from matplotlib<4.0.0,>=3.6.1->epimodels==0.4.0) (0.11.0)
Requirement already satisfied: pillow>=6.2.0 in /usr/lib/python3/dist-packages (from matplotlib<4.0.0,>=3.6.1->epimodels==0.4.0) (9.0.1)
Requirement already satisfied: fonttools>=4.22.0 in /usr/lib/python3/dist-packages (from matplotlib<4.0.0,>=3.6.1->epimodels==0.4.0) (4.29.1)
Collecting contourpy>=1.0.1
  Using cached contourpy-1.0.5-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (295 kB)
Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.10/dist-packages (from matplotlib<4.0.0,>=3.6.1->epimodels==0.4.0) (2.8.2)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/lib/python3/dist-packages (from matplotlib<4.0.0,>=3.6.1->epimodels==0.4.0) (1.3.2)
Requirement already satisfied: pyparsing>=2.2.1 in /usr/lib/python3/dist-packages (from matplotlib<4.0.0,>=3.6.1->epimodels==0.4.0) (2.4.7)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.1.2->salib) (2022.2.1)
Requirement already satisfied: mpmath>=0.19 in /usr/local/lib/python3.10/dist-packages (from sympy<2.0.0,>=1.11.1->epimodels==0.4.0) (1.2.1)
Requirement already satisfied: dill>=0.3.5.1 in /usr/local/lib/python3.10/dist-packages (from multiprocess->salib) (0.3.5.1)
Requirement already satisfied: scikit-learn>=0.16.0 in /usr/local/lib/python3.10/dist-packages (from pyitlib->epimodels==0.4.0) (1.1.1)
Requirement already satisfied: future>=0.16.0 in /usr/lib/python3/dist-packages (from pyitlib->epimodels==0.4.0) (0.18.2)
Requirement already satisfied: six>=1.5 in /usr/lib/python3/dist-packages (from python-dateutil>=2.7->matplotlib<4.0.0,>=3.6.1->epimodels==0.4.0) (1.16.0)
Requirement already satisfied: joblib>=1.0.0 in /usr/local/lib/python3.10/dist-packages (from scikit-learn>=0.16.0->pyitlib->epimodels==0.4.0) (1.1.0)
Requirement already satisfied: threadpoolctl>=2.0.0 in /usr/local/lib/python3.10/dist-packages (from scikit-learn>=0.16.0->pyitlib->epimodels==0.4.0) (3.1.0)
Installing collected packages: sympy, numpy, mypy, cython, scipy, contourpy, matplotlib, pyitlib, epimodels
  WARNING: The script isympy is installed in '/home/fccoelho/.sage/local/bin' which is not on PATH.
  Consider adding this directory to PATH or, if you prefer to suppress this warning, use --no-warn-script-location.
  WARNING: The scripts f2py, f2py3 and f2py3.10 are installed in '/home/fccoelho/.sage/local/bin' which is not on PATH.
  Consider adding this directory to PATH or, if you prefer to suppress this warning, use --no-warn-script-location.
  WARNING: The scripts dmypy, mypy, mypyc, stubgen and stubtest are installed in '/home/fccoelho/.sage/local/bin' which is not on PATH.
  Consider adding this directory to PATH or, if you prefer to suppress this warning, use --no-warn-script-location.
  WARNING: The scripts cygdb, cython and cythonize are installed in '/home/fccoelho/.sage/local/bin' which is not on PATH.
  Consider adding this directory to PATH or, if you prefer to suppress this warning, use --no-warn-script-location.
  Attempting uninstall: epimodels
    Found existing installation: epimodels 0.3.21
    Uninstalling epimodels-0.3.21:
      Successfully uninstalled epimodels-0.3.21
ERROR: pip's dependency resolver does not currently take into account all the packages that are installed. This behaviour is the source of the following dependency conflicts.
pysus 0.6.0 requires numpy==1.23.2, but you have numpy 1.23.3 which is incompatible.
numba 0.56.0 requires numpy<1.23,>=1.18, but you have numpy 1.23.3 which is incompatible.
Successfully installed contourpy-1.0.5 cython-0.29.32 epimodels-0.4.0 matplotlib-3.6.1 mypy-0.982 numpy-1.23.3 pyitlib-0.2.2 scipy-1.9.2 sympy-1.11.1

%display typeset

from SALib.sample import saltelli
from SALib.analyze import sobol
from epimodels.continuous import models as cm
import seaborn as sns
import numpy as np
import pandas as pd
import pylab as plt

Escolhendo o Modelo

Vamos usar o Modelo SEIAHR apresentado no suplemento 2 e disponível na biblioteca epimodels.

my_model = cm.SEQIAHR()

my_model.parameters

$$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|OrderedDict([('chi',|\verb| |\verb|'$\chiParseError: KaTeX parse error: \verb ended by end of line instead of matching delimiter\phiParseError: KaTeX parse error: \verb ended by end of line instead of matching delimiter\betaParseError: KaTeX parse error: \verb ended by end of line instead of matching delimiter\rhoParseError: KaTeX parse error: \verb ended by end of line instead of matching delimiter\deltaParseError: KaTeX parse error: \verb ended by end of line instead of matching delimiter\gammaParseError: KaTeX parse error: \verb ended by end of line instead of matching delimiter\alphaParseError: KaTeX parse error: \verb ended by end of line instead of matching delimiter\muParseError: KaTeX parse error: \verb ended by end of line instead of matching delimiterpParseError: KaTeX parse error: \verb ended by end of line instead of matching delimiterqParseError: KaTeX parse error: \verb ended by end of line instead of matching delimiterr$')])|$$

Vamos atribuir valores iniciais aos parâmetros do modelo

params = {
'chi': .3,
'phi': 0.012413633926076584,
'beta': 0.27272459855759813,
'rho': 0.2190519831830368,
'delta': 0.04168480042146949,
'gamma': 0.04,
'alpha':  0.3413355572047603,
'mu': 0.02359234606623134,
'p': 0.7693029079871165,
'q': 50,
'r': 55,
}
inits = [.99, 0, 1e-4, 0, 0, 0, 0, 0]

Durante a análise de sensibilidade vai ser necessário simular muitas vezes o modelo, então é interessante termos uma ideia de quanto tempo custa simulá-lo

%%timeit
# parms = dict(zip(params.keys(),param_values[3]))
my_model(inits=inits, trange = [0,200], totpop=1, params=params)
# max(my_model.traces['I'])

32.4 ms ± 591 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

my_model.plot_traces(['I','A'])

Definindo a análise

problem = {
    'num_vars': 11,
    'names': list(params.keys()),
    'bounds': [[0,.3],[0,1.],[0.1,2],[0.2,1],[0.01,1],[0.01,1],
              [0.01,1],[0.01,.7],[0.01,.7],[1,30],[1,12]]
}

Gerando as amostras

param_values = saltelli.sample(problem, 1024)

/tmp/ipykernel_3634413/1862108280.py:1: DeprecationWarning: `salib.sample.saltelli` will be removed in SALib 1.5. Please use `salib.sample.sobol`
  param_values = saltelli.sample(problem, Integer(1024))

param_values.shape

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(24576, 11\right)$$

(32.4*24576)/1000

$$\renewcommand{\Bold}[1]{\mathbf{#1}}796.262400000000$$

Rodando o Modelo

Para fazer a análise de sensibilidade precisamos selecionar um aspecto da saída do model, sobre o qual desejamos estudar a variância em resposta à varância dos parâmetros. Para este caso simples, vamos escolher o valor de pico de $I(t)$.

def eval_model(parametros):
    parms = dict(zip(params.keys(),parametros))
    mod = cm.SEQIAHR()
    mod(inits=inits, trange = [0,200], totpop=1, params=parms)
    Y = max(mod.traces['I'])
    return Y

Agora, basta executar uma simulação para cada conjunto de valores de parâmetros que amostramos.

# Y = np.zeros([param_values.shape[0]])
# for i, p in enumerate(param_values):
#     Y[i] = eval_model(p)

Para acelerar a execução dos cenários, vamos paralelizar a simulação do modelo

from multiprocessing import Pool

Po = Pool()
Y = Po.map(eval_model, param_values)

Po.close()

Si = sobol.analyze(problem, np.array(Y), print_to_console=False)

def plot_sobol(si,prob, order=1):
    Si_filter = {k:si[k] for k in ['ST','ST_conf','S1','S1_conf']}
    Si_df = pd.DataFrame(Si_filter, index=problem['names'][:-1])
    fig, ax = plt.subplots(1, figsize=(10,8))

    indices = Si_df[['S1','ST']]
    err = Si_df[['S1_conf','ST_conf']]

    indices.plot.bar(yerr=err.values.T,ax=ax)
plot_sobol(Si,problem)

Para plotar os índices de segunda ordem, vamos usar uma outra abordagem para a visualização

import itertools
from math import pi
from matplotlib.legend_handler import HandlerPatch


def normalize(x, xmin, xmax):
    return (x-xmin)/(xmax-xmin)


def plot_circles(ax, locs, names, max_s, stats, smax, smin, fc, ec, lw,
                 zorder):
    s = np.asarray([stats[name] for name in names])
    s = 0.01 + max_s * np.sqrt(normalize(s, smin, smax))

    fill = True
    for loc, name, si in zip(locs, names, s):
        if fc=='w':
            fill=False
        else:
            ec='none'

        x = np.cos(loc)
        y = np.sin(loc)

        circle = plt.Circle((x,y), radius=si, ec=ec, fc=fc, transform=ax.transData._b,
                            zorder=zorder, lw=lw, fill=True)
        ax.add_artist(circle)


def filter(sobol_indices, names, locs, criterion, threshold):
    if criterion in ['ST', 'S1', 'S2']:
        data = sobol_indices[criterion]
        data = np.abs(data)
        data = data.flatten() # flatten in case of S2
        # TODO:: remove nans

        filtered = ([(name, locs[i]) for i, name in enumerate(names[:-1]) if
                     data[i]>threshold])
        filtered_names, filtered_locs = zip(*filtered)
    elif criterion in ['ST_conf', 'S1_conf', 'S2_conf']:
        raise NotImplementedError
    else:
        raise ValueError('unknown value for criterion')

    return filtered_names, filtered_locs


def plot_sobol_indices(sobol_indices, criterion='ST', threshold=0.01):
    '''plot sobol indices on a radial plot

    Parameters
    ----------
    sobol_indices : dict
                    the return from SAlib
    criterion : {'ST', 'S1', 'S2', 'ST_conf', 'S1_conf', 'S2_conf'}, optional
    threshold : float
                only visualize variables with criterion larger than cutoff

    '''
    max_linewidth_s2 = 15#25*1.8
    max_s_radius = 0.3

    # prepare data
    # use the absolute values of all the indices
    #sobol_indices = {key:np.abs(stats) for key, stats in sobol_indices.items()}

    # dataframe with ST and S1
    sobol_stats = {key:sobol_indices[key] for key in ['ST', 'S1']}
    sobol_stats = pd.DataFrame(sobol_stats, index=problem['names'][:-1])

    smax = sobol_stats.max().max()
    smin = sobol_stats.min().min()

    # dataframe with s2
    s2 = pd.DataFrame(sobol_indices['S2'], index=problem['names'][:-1],
                      columns=problem['names'][:-1])
    s2[s2<0.0]=0. #Set negative values to 0 (artifact from small sample sizes)
    s2max = s2.max().max()
    s2min = s2.min().min()

    names = problem['names'][:-1]
    n = len(names)
    ticklocs = np.linspace(0, 2*pi, n)
    locs = ticklocs[0:-1]

    filtered_names, filtered_locs = filter(sobol_indices, names, locs,
                                           criterion, threshold)

    # setup figure
    fig = plt.figure(figsize=(15,15))
    ax = fig.add_subplot(111, polar=True)
    ax.grid(False)
    ax.spines['polar'].set_visible(False)
    ax.set_xticks(ticklocs)

    ax.set_xticklabels(names, fontsize=24)
    ax.set_yticklabels([])
    ax.set_ylim(top=1.4)
    legend(ax)

    # plot ST
    plot_circles(ax, filtered_locs, filtered_names, max_s_radius,
                 sobol_stats['ST'], smax, smin, 'w', 'k', 1, 9)

    # plot S1
    plot_circles(ax, filtered_locs, filtered_names, max_s_radius,
                 sobol_stats['S1'], smax, smin, 'k', 'k', 1, 10)

    # plot S2
    for name1, name2 in itertools.combinations(zip(filtered_names, filtered_locs), 2):
        name1, loc1 = name1
        name2, loc2 = name2

        weight = s2.loc[name1, name2]
        lw = 0.5+max_linewidth_s2*normalize(weight, s2min, s2max)
        ax.plot([loc1, loc2], [1,1], c='darkgray', lw=lw, zorder=1)

    return fig



class HandlerCircle(HandlerPatch):
    def create_artists(self, legend, orig_handle,
                       xdescent, ydescent, width, height, fontsize, trans):
        center = 0.5 * width - 0.5 * xdescent, 0.5 * height - 0.5 * ydescent
        p = plt.Circle(xy=center, radius=orig_handle.radius)
        self.update_prop(p, orig_handle, legend)
        p.set_transform(trans)
        return [p]

def legend(ax):
    some_identifiers = [plt.Circle((0,0), radius=5, color='k', fill=False, lw=1),
                        plt.Circle((0,0), radius=5, color='k', fill=True),
                        plt.Line2D([0,0.5], [0,0.5], lw=8, color='darkgray')]
    ax.legend(some_identifiers, ['ST', 'S1', 'S2'],
              loc=(1,0.75), borderaxespad=0.1, mode='expand', fontsize=16,
              handler_map={plt.Circle: HandlerCircle()})


sns.set_style('whitegrid')
fig = plot_sobol_indices(Si, criterion='ST', threshold=0.005)