stats

QSDsan has many functions uncertainty analysis visualization and sensitivity analysis/visualization. You can find the complete script of following code snippets from stats_demo.py in the bwaise module of Exposan.

Uncertainties

stats.plot_uncertainties(x_axis=(), y_axis=(), kind='box', adjust_hue=False, file='', close_fig=True, center_kws={}, margin_kws={})

Visualize uncertainty analysis results as one of the following depending on inputs:

input

returned plot

x_axis

y_axis

kind

orientation

center

margin

single or iterable

None

box

1D horizontal

box

N/A

hist

histogram

kde

kernel density

None

single or iterable

box

1D vertical

box

N/A

hist

histogram

kde

kernel density

single

single

hist-box

2D

histogram

box

hist-kde

kernel density

hist-hist

histogram

kde-box

kernel density

box

kde-hist

histogram

kde-kde

kernel density

Note

When both x_axis and y_axis are not None (i.e., the figure is 2D), for clarity reasons, this function requires len(x_axis) and len(y_axis) to be both 1.

If wanted, it is possible to use colors (hue in seaborn) to differentiated the different parameters or metrics in either x_axis or y_axis. To achieve this, please refer to the documentation of seaborn.

Similarly, there are other potential combinations of the center and margin plots that are supported by seaborn but not included here.

Parameters:

  • model (biosteam.Model) – The model with uncertainty analysis (in <Model.table>) results for plotting.

  • x_axis (biosteam.Parameter, biosteam.Metric or Iterable) – What to plot on the x-axis, can be parameters or metrics of the model, default to all model metrics included in the model result table if neither x nor y is provided.

  • y_axis (biosteam.Parameter, biosteam.Metric or Iterable) – What to plot on the y-axis, can be parameters or metrics of the model, default to None.

  • kind (str) – What kind of plot to be returned, refer to the summary table for valid inputs.

  • adjust_hue (bool) – Whether to adjust the hue of the colors based on the data.

  • file (str) – If provided, the generated figure will be saved as a png file.

  • close_fig (bool) – Whether to close the figure (if not close, new figure will be overlaid on the current figure).

  • center_kws (dict) – Will be passed to seaborn for the center plot.

  • margin_kws (dict) – Will be passed to seaborn for the margin plots.

Returns:

  • figure (matplotlib.figure.Figure) – The generated figure.

  • axis (matplotlib.axes._subplots.AxesSubplot or Iterable) – The generated figure axis (or axes for 2D figure).

Examples

qsdsan.stats

See also

get_correlations(), seaborn.jointplot()

Examples

Box plot

import pandas as pd
from qsdsan import stats as s
from exposan import bwaise as bw

m = bw.models
modelA = bw.modelA

# Total COD/N/P/K recovery and net cost/GWP
modelA.metrics = key_metrics = bw.get_key_metrics(
    modelA, alt_names={'Annual net cost': 'Cost',
                       'Net emission GlobalWarming': 'GWP'})

seed = 3221 # set numpy seed for sample reproducibility

# Run Monte Carlo uncertainty analysis and get Spearman rank correlations,
# here we use a small sample size for demonstrative purpose
m.run_uncertainty(modelA, N=100, seed=seed, rule='L',
                  percentiles=(0, 0.05, 0.25, 0.5, 0.75, 0.95, 1))

# Pass a path to `file` or use `fig.savefig` if want to save the figure,
# the `file` kwarg exists for pretty much all of the plotting functions
fig, ax = s.plot_uncertainties(modelA,
                               x_axis=key_metrics[:-2], # only recoveries
                               kind='box', file='')

# Trim figure
fig.subplots_adjust(bottom=0.25)
for label in ax.get_xticklabels():
    label.set_rotation(45)
../_images/plot_uncer_box.png

Histogram plot

# Kernel density curve can be added to the histogram,
# with a log scale, we can have all metric results in the same plot
fig, ax = s.plot_uncertainties(modelA, y_axis=key_metrics, kind='hist',
                               center_kws={'kde':True, 'log_scale': 10})
../_images/plot_uncer_hist.png 
# We can also have 2D histogram plot
fig, axes = s.plot_uncertainties(modelA,
                                 x_axis=key_metrics[-2], # cost
                                 y_axis=key_metrics[-1], # GWP
                                 kind='hist-box')
../_images/plot_uncer_hist-box.png

Kernel density plots

# Similar to histogram plots, kernel density plots can be 1D
fig, ax = s.plot_uncertainties(modelA, x_axis=key_metrics, kind='kde',
                               center_kws={'fill': True, 'log_scale': 2})
../_images/plot_uncer_kde.png 
# Or 2D with different kinds of margins
fig, axes = s.plot_uncertainties(modelA, x_axis=key_metrics[-2],
                                 y_axis=key_metrics[-1], kind='kde-kde',
                                 center_kws={'fill': True})
../_images/plot_uncer_kde-kde.png 
fig, axes = s.plot_uncertainties(modelA, x_axis=key_metrics[-2],
                                 y_axis=key_metrics[-1], kind='kde-hist',
                                 center_kws={'fill': True},
                                 margin_kws={'kde': True, 'fill': False})
../_images/plot_uncer_kde-hist.png

Correlations

stats.get_correlations(input_x=None, input_y=None, kind='Spearman', nan_policy='propagate', file='', **kwargs)

Get correlation coefficients between two inputs using scipy.

Parameters:

  • model (biosteam.Model) – Uncertainty model with defined parameters and metrics.

  • input_x (biosteam.Parameter or biosteam.Metric) – The first set of input, can be single values or an iterable, will be defaulted to all model parameters if not provided.

  • input_y (biosteam.Parameter or biosteam.Metric) – The second set of input, can be single values or an iterable, will be defaulted to all model parameters if not provided.

  • kind (str) – Can be “Spearman” for Spearman’s rho, “Pearson” for Pearson’s r, “Kendall” for Kendall’s tau, or “KS” for Kolmogorov–Smirnov’s D.

  • nan_policy (str) –

    • “propagate”: returns nan.

    • ”raise”: raise an error.

    • ”omit”: drop the pair from analysis.

  • file (str) – If provided, the results will be saved as an Excel file.

  • **kwargs (dict) – Other kwargs that will be passed to scipy.

Return type:

Two pandas.DataFrame containing the test correlations and p-values.

Examples

Run uncertainty and sensitivity analyses

>>> # Import a pre-constructed system model
>>> from qsdsan import stats as s
>>> from qsdsan.utils import create_example_model
>>> # Use a small sample size for demonstrative purpose
>>> model = create_example_model(evaluate=True, N=100, rule='L', seed=554)
>>> p, m = model.parameters, model.metrics
>>> # Correlations, p-values
>>> r_df1, p_df1 = s.get_correlations(model, kind='Spearman')
>>> r_df2, p_df2 = s.get_correlations(model, kind='Pearson')
>>> r_df3, p_df3 = s.get_correlations(model, kind='Kendall')
>>> r_df4, p_df4 = s.get_correlations(
...     model, input_x=p[0], input_y=m[-1], kind='KS', thresholds=[0.1])

Plots for uncertainty analysis

>>> # Box
>>> fig, ax = s.plot_uncertainties(model, x_axis=m[0], kind='box')
>>> # Histogram
>>> fig, ax = s.plot_uncertainties(model, y_axis=m[1], kind='hist')
>>> # Kernel density
>>> fig, ax = s.plot_uncertainties(model, x_axis=m[2], kind='kde')
>>> # 2-D counterparts
>>> fig, ax = s.plot_uncertainties(
...     model, x_axis=m[0], y_axis=m[1], kind='hist-hist')
>>> fig, ax = s.plot_uncertainties(
...     model, x_axis=m[0], y_axis=m[1], kind='kde-kde')

Plots for sensitivity analysis

>>> fig, ax = s.plot_correlations(r_df1, metrics=m[-2])
>>> fig, ax = s.plot_correlations(r_df2)

See also

qsdsan.stats

scipy.stats.spearmanr()

scipy.stats.pearsonr()

scipy.stats.kendalltau()

scipy.stats.kstest()

stats.plot_correlations(parameters=(), metrics=(), top=None, file='', close_fig=True, **kwargs)

Visualize the correlations between model parameters and metric results as tornado (single metric) or bubble plots (multiple metrics).

Parameters:

  • result_df (pandas.DataFrame) – Result table generated by get_correlations() containing correlation indices.

  • parameters (biosteam.Parameter) – Metric(s) of interest for the plot, will be default to all parameters included in result_df if not provided.

  • metrics (biosteam.Metric) – Metric(s) of interest for the plot, will be default to all metrics included in result_df if not provided.

  • top (int) – Plot the top X parameters with the highest absolute correlation indices, this is only applicable for the case of just one metric.

  • file (str) – If provided, the generated figure will be saved as a png file.

  • close_fig (bool) – Whether to close the figure (if not close, new figure will be overlaid on the current figure).

  • **kwargs (dict) – Other kwargs that will be passed to seaborn.relplot().

Returns:

  • figure (matplotlib.figure.Figure) – The generated figure.

  • axis (matplotlib.axes._subplots.AxesSubplot) – The generated figure axis.

Examples

get_correlations()

qsdsan.stats

See also

biosteam.plots.plot_spearman(), seaborn.relplot()

Examples

Bar plot for single metric

spearman_rho, spearman_p = s.get_correlations(
    modelA, kind='Spearman', nan_policy='raise',
    file='') # pass a path to `file` if you want to save the results as an Excel

# Filter out parameters that only meet a certain threshold
def filter_parameters(model, df, threshold):
    new_df = pd.concat((df[df>=threshold], df[df<=-threshold]))
    filtered = new_df.dropna(how='all')
    param_dct = {p.name_with_units:p for p in model.get_parameters()}
    parameters = set(param_dct[i[1]] for i in filtered.index)
    return list(parameters)

# Only want parameters with Spearman's rho >= 0.4 or <= -0.4
modelA.parameters = key_parameters = \
    filter_parameters(modelA, spearman_rho, threshold=0.4)

fig, ax = s.plot_correlations(spearman_rho, parameters=key_parameters,
                                  metrics=key_metrics[-2])

fig.subplots_adjust(left=0.25)
../_images/plot_corr_bar.png

Bubble plot for multiple metrics

fig, ax = s.plot_correlations(
    spearman_rho, parameters=key_parameters, metrics=key_metrics)
../_images/plot_corr_bubble.png

Input and sample preparation

stats.define_inputs()

Define the model inputs (referred to as “problem”) to be used for sampling by SALib.

Parameters:

model (biosteam.Model) – Uncertainty model with defined parameters and metrics.

Returns:

inputs – A dict containing model inputs for sampling by SALib.

Return type:

dict

See also

qsdsan.stats

SALib Basics

stats.generate_samples(kind, N, seed=None, **kwargs)

Generate samples for sensitivity analysis using SALib.

Parameters:

  • model (biosteam.Model) – Uncertainty model with defined parameters and metrics.

  • inputs (dict) – A dict generated by qsdsan.sensitivity.define_inputs() to be used for SALib, keys should include “num_vars”, “names”, and “bounds”.

  • kind (str) – Can be “FAST” for Fourier amplitude sensitivity test (FAST), “RBD” for random balance design FAST (latin hypercube sampling [LHS]), “Morris” for Morris One-at-A-Time (OAT), or “Sobol” for Sobol analysis (Saltelli sampling).

  • N (int) – The number of samples or trajectories (Morris).

  • seed (int) – Seed to generate random samples.

  • **kwargs (dict) – Other kwargs that will be passed to SALib.

Returns:

samples – Samples to be used for the indicated sensitivity analyses.

Return type:

array

See also

qsdsan.stats

SALib API

Morris

stats.morris_analysis(inputs, metrics=None, nan_policy='propagate', conf_level=0.95, print_to_console=False, file='', **kwargs)

Run Morris sensitivity analysis using SALib.

Parameters:

  • model (biosteam.Model) – Uncertainty model with defined parameters and metrics.

  • inputs (dict) – A dict generated by qsdsan.sensitivity.define_inputs() to be used for SALib, keys should include “num_vars”, “names”, and “bounds”.

  • metrics (biosteam.Metric) – Metrics to be included for Morris analysis, must be a subset of the metrics of the model to be analyzed. (i.e., included in the metrics attribute of the given model). If None is provided, all metrics in the model will be included.

  • nan_policy (str) –

    • “propagate”: returns nan.

    • ”raise”: raise an error.

    • ”fill_mean”: fill nan with mean of the results.

  • conf_level (float) – Confidence level of results.

  • print_to_console (bool) – Whether to show results in the console.

  • file (str) – If provided, the results will be saved as an Excel file.

  • **kwargs (dict) – Other kwargs that will be passed to SALib.

Returns:

morris_dct – A dict of Morris analysis results.

Return type:

dict

Examples

>>> # Import a pre-constructed system model
>>> from qsdsan import stats as s
>>> from qsdsan.utils import create_example_model
>>> model = create_example_model(evaluate=False)
>>> # Morris analysis requires special samples
>>> inputs = s.define_inputs(model)
>>> # Use a small sample size for demonstrative purpose
>>> samples = s.generate_samples(inputs, kind='Morris', N=10, seed=554)
>>> model.load_samples(samples)
>>> model.evaluate()
>>> dct = s.morris_analysis(model, inputs, seed=554, nan_policy='fill_mean')
>>> fig, ax = s.plot_morris_results(dct, metric=model.metrics[0])

See also

qsdsan.stats

SALib API

stats.morris_till_convergence(inputs, metrics=None, N_max=20, seed=None, threshold=0.1, nan_policy='propagate', conf_level=0.95, print_to_console=False, file='', **kwargs)

Run Morris analysis from N=2 to N=N_max until the results converge (i.e., mu_star_conf/mu_star_max < threshold for all parameters, where as mu_star_max is the maximum \({\mu^*}\) value for a certain metric, and this should be satisfied for all metrics).

Parameters:

  • model (biosteam.Model) – Uncertainty model with defined parameters and metrics.

  • inputs (dict) – A dict generated by qsdsan.sensitivity.define_inputs() to be used for SALib, keys should include “num_vars”, “names”, and “bounds”.

  • metrics (biosteam.Metric) – Metrics to be included for Morris analysis, must be a subset of the metrics of the model to be analyzed. (i.e., included in the metrics attribute of the given model). If None is provided, all metrics in the model will be included.

  • N_max (int) – Maximum number of trajectories to be considered.

  • seed (int) – Seed to generate random samples.

  • threshold (float) – Threshold for the convergence.

  • nan_policy (str) –

    • “propagate”: returns nan.

    • ”raise”: raise an error.

    • ”fill_mean”: fill nan with mean of the results.

  • conf_level (float) – Confidence level of results.

  • print_to_console (bool) – Whether to show results in the console.

  • file (str) – If provided, the results will be saved as an Excel file.

  • **kwargs (dict) – Other kwargs that will be passed to SALib.

Examples

>>> # Import a pre-constructed system model
>>> from qsdsan import stats as s
>>> from qsdsan.utils import create_example_model
>>> model = create_example_model(evaluate=False)
>>> # Morris analysis requires special samples
>>> inputs = s.define_inputs(model)
>>> # Use a small maximum trajectory number for demonstrative purpose
>>> dct = s.morris_till_convergence(model, inputs, seed=554, N_max=10)
mu_star has not converged within 10 trajectories.
>>> fig, ax = s.plot_morris_convergence(dct, metric=model.metrics[-2], plot_rank=False)
>>> fig, ax = s.plot_morris_convergence(dct, metric=model.metrics[-2], plot_rank=True)

See also

qsdsan.stats

qsdsan.stats.morris_analysis()

stats.plot_morris_results(metric, kind='scatter', ax=None, x_axis='mu_star', plot_lines=True, k1=0.1, k2=0.5, k3=1, label_kind='number', color='k', file='', close_fig=True, **kwargs)

Visualize the results from Morris One-at-A-Time analysis as either scatter or bar plots. In scatter plots, the x values are \({\mu^*}\) and the y values are \({\sigma}\). In bar plots, bar length indicate the \({\mu^*}\) values with error bars representing confidence intervals of the analysis.

Parameters:

  • morris_dct (dict) – Results dict generated by morris_analysis().

  • metric (biosteam.Metric) – The metric of interest for the plot.

  • kind (str) – Either “scatter” (\({\sigma}\)) vs. \({\mu^*}\) or “bar” (\({\mu^*}\) with confidence interval) plot.

  • ax (matplotlib.AxesSubplot) – If provided, the figure will be plotted on this axis.

  • x_axis (str) – X-axis parameter, should be either “mu_star” (the commonly used one) or “mu”.

  • k1 (float) – The slope to differentiate monotonic (above the line) and linear (below the line).

  • k2 (float) – The slope to differentiate almost monotonic (above the line) and monotonic (below the line).

  • k3 (float) – The slope to differentiate non-linear and/or non-monotonic (above the line) and almost monotonic (below the line).

  • label_kind (str) – How to label the points in the scatter plot, can be “number” (use index number of the result table), “name” (use index name of the result table), or None (no labels).

  • color (str or RGBs) – Plot color.

  • file (str) – If provided, the generated figure will be saved as a png file.

  • close_fig (bool) – Whether to close the figure (if not close, new figure will be overlaid on the current figure).

  • **kwargs (dict) – Other kwargs that will be passed to morris.horizontal_bar_plot() in SALib.plotting.

Returns:

  • figure (matplotlib.figure.Figure) – The generated figure.

  • axis (matplotlib.axes._subplots.AxesSubplot) – The generated figure axis.

See also

morris_analysis()

qsdsan.stats

stats.plot_morris_convergence(metric, parameters=(), plot_rank=False, kind='line', ax=None, show_error=True, palette='pastel', file='', close_fig=True)

Plot the evolution of \({\mu^*}\) or its rank with the number of trajectories.

Parameters:

  • result_dct (dict) – Result dictionary generated from qsdsan.stats.morris_till_convergence()

  • metric (biosteam.Metric) – The metric of interest for the plot.

  • parameters (biosteam.Parameter) – Single or a Iterable of model parameters whose \({\mu^*}\) will be included in the plot. Will be set to all parameters in retult_dct will be used if not provided.

  • plot_rank (bool) –

    If True, will plot rank of \({\mu^*}\) instead of value.

    Note

    If plot_rank is True, error bars or bands will not be included.

  • kind (str) – Can be either ‘line’ or ‘scatter’.

  • ax (matplotlib.AxesSubplot) – If provided, the figure will be plotted on this axis.

  • show_error (bool) – Whether to include the confidence interval in the plot, will be bars for scatter plot and bands for line plot.

  • palette (string, list, dict, or matplotlib.colors.Colormap) – Will be passed on to seaborn.color_palette().

  • file (str) – If provided, the generated figure will be saved as a png file.

  • close_fig (bool) – Whether to close the figure (if not close, new figure will be overlaid on the current figure).

Returns:

  • figure (matplotlib.figure.Figure) – The generated figure.

  • axis (matplotlib.axes._subplots.AxesSubplot) – The generated figure axis.

See also

morris_analysis()

qsdsan.stats

Examples

\(\sigma\) vs. \(\mu^*\)

# Run Morris analysis without testing the convergence,
# here we use a small sample size for demonstrative purpose
inputs = s.define_inputs(modelA)
morris_samples = s.generate_samples(inputs, kind='Morris', N=10, seed=seed)

evaluate = bw.evaluate
evaluate(modelA, morris_samples)

dct = s.morris_analysis(modelA, inputs, metrics=key_metrics, seed=seed,
                        nan_policy='fill_mean')

# Unfortunately the auto-labelling is not good when you have close points,
# so you'll have to do some manual manipulation
fig, ax = s.plot_morris_results(dct, key_metrics[-2])

fig.subplots_adjust(bottom=0.3)
../_images/plot_morris.png

Line plot with error bands for evolutionary of \(\mu^*\)

# Test if mu_star can converge within 100 trajectories
# (spoiler: it cannot because we already sort of selected the key parameters,
# and you will get a message prompt)
dct = s.morris_till_convergence(modelA, inputs, metrics=key_metrics, seed=seed,
                                N_max=100)

# Look at mu_star values for two parameters with regard to cost
fig, ax = s.plot_morris_convergence(dct,
                                    parameters=key_parameters[:2],
                                    metric=key_metrics[-2], plot_rank=False)
../_images/plot_morris_conv.png

Line plot for evolutionary of \(\mu^*\) rank

# Look at ranks of mu_star values for all parameters with regard to cost
fig, ax = s.plot_morris_convergence(dct, parameters=key_parameters,
                                    metric=key_metrics[-2], plot_rank=True)
../_images/plot_morris_conv_rank.png

FAST

stats.fast_analysis(inputs, kind, metrics=None, nan_policy='propagate', conf_level=0.95, print_to_console=False, file='', **kwargs)

Run Fourier amplitude sensitivity test (Saltelli’s extended FAST) or random balance design (RBD) FAST using SALib.

Parameters:

  • model (biosteam.Model) – Uncertainty model with defined parameters and metrics.

  • inputs (dict) – A dict generated by qsdsan.sensitivity.define_inputs() to be used for SALib, keys should include “num_vars”, “names”, and “bounds”.

  • kind (str) – Either “FAST” or “RBD”.

  • metrics (biosteam.Metric) – Metrics to be included for Morris analysis, must be a subset of the metrics of the model to be analyzed. (i.e., included in the metrics attribute of the given model). If None is provided, all metrics in the model will be included.

  • nan_policy (str) –

    • “propagate”: returns nan.

    • ”raise”: raise an error.

    • ”fill_mean”: fill nan with mean of the results.

  • conf_level (float) – Confidence level of results.

  • print_to_console (bool) – Whether to show results in the console.

  • file (str) – If provided, the results will be saved as an Excel file.

  • **kwargs (dict) – Other kwargs that will be passed to SALib.

Returns:

fast_dct – A dict of FAST analysis results.

Return type:

dict

Examples

>>> # Import a pre-constructed system model
>>> from qsdsan import stats as s
>>> from qsdsan.utils import create_example_model
>>> model = create_example_model(evaluate=False)
>>> # FAST analysis requires special samples
>>> inputs = s.define_inputs(model)
>>> # Use a small sample size for demonstrative purpose
>>> samples = s.generate_samples(inputs, kind='FAST', N=100, M=4, seed=554)
>>> model.load_samples(samples)
>>> model.evaluate()
>>> dct = s.fast_analysis(model, inputs, kind='FAST', M=4, seed=554, nan_policy='fill_mean')
>>> fig, ax = s.plot_fast_results(dct, metric=model.metrics[-3])
>>> # If want to do RBD-FAST
>>> # Use a small sample size for demonstrative purpose
>>> samples = s.generate_samples(inputs, kind='RBD', N=100, seed=554)
>>> model.load_samples(samples)
>>> model.evaluate()
>>> dct = s.fast_analysis(model, inputs, kind='RBD', seed=554, nan_policy='fill_mean')
>>> fig, ax = s.plot_fast_results(dct, metric=model.metrics[-3])

See also

qsdsan.stats

SALib API

stats.plot_fast_results(metric, parameters=(), ax=None, error_bar=True, file='', close_fig=True)

Visualize the results from FAST or RBD-FAST analysis as a bar plot.

Parameters:

  • result_dct (dict) – Result dictionary generated from qsdsan.stats.fast_analysis()

  • metric (biosteam.Metric) – The metric of interest for the plot.

  • parameters (biosteam.Parameter) – Single or a Iterable of model parameters whose \({\mu^*}\) will be included in the plot. Will be set to all parameters in retult_dct will be used if not provided.

  • ax (matplotlib.AxesSubplot) – If provided, the figure will be plotted on this axis.

  • error_bar (bool) – Whether to include the confidence interval as error bars in the plot.

  • file (str) – If provided, the generated figure will be saved as a png file.

  • close_fig (bool) – Whether to close the figure (if not close, new figure will be overlaid on the current figure).

Returns:

  • figure (matplotlib.figure.Figure) – The generated figure.

  • axis (matplotlib.axes._subplots.AxesSubplot) – The generated figure axis.

See also

fast_analysis()

qsdsan.stats

Examples

Bar plot for (e)FAST

# Total and main effects from FAST analysis,
# here we use a small sample size for demonstrative purpose
fast_samples = s.generate_samples(inputs, kind='FAST', N=100, M=4, seed=seed)

evaluate(modelA, fast_samples)

dct = s.fast_analysis(modelA, inputs, kind='FAST', metrics=key_metrics,
                      M=4, seed=seed, nan_policy='fill_mean')

fig, ax = s.plot_fast_results(dct, metric=key_metrics[-2])

fig.subplots_adjust(left=0.4)
../_images/plot_fast.png

Bar plot for RBD-FAST

# Main effects from RBD-FAST analysis,
# here we use a small sample size for demonstrative purpose
fast_samples = s.generate_samples(inputs, kind='RBD', N=100, seed=seed)

evaluate(modelA, fast_samples)

dct = s.fast_analysis(modelA, inputs, kind='RBD', metrics=key_metrics,
                      seed=seed, nan_policy='fill_mean')

fig, ax = s.plot_fast_results(dct, metric=key_metrics[-2])

fig.subplots_adjust(left=0.4)
../_images/plot_rbd.png

Sobol

stats.sobol_analysis(inputs, metrics=None, nan_policy='propagate', calc_second_order=True, conf_level=0.95, print_to_console=False, file='', **kwargs)

Run Sobol sensitivity analysis using SALib.

Parameters:

  • model (biosteam.Model) – Uncertainty model with defined parameters and metrics.

  • inputs (dict) – A dict generated by qsdsan.sensitivity.define_inputs() to be used for SALib, keys should include “num_vars”, “names”, and “bounds”.

  • metrics (biosteam.Metric) – Metrics to be included for Morris analysis, must be a subset of the metrics of the model to be analyzed. (i.e., included in the metrics attribute of the given model). If None is provided, all metrics in the model will be included.

  • nan_policy (str) –

    • “propagate”: returns nan.

    • ”raise”: raise an error.

    • ”fill_mean”: fill nan with mean of the results.

  • calc_second_order (bool) – Whether to calculate second-order interaction effects.

  • conf_level (float) – Confidence level of results.

  • print_to_console (bool) – Whether to show results in the console.

  • file (str) – If provided, the results will be saved as an Excel file.

  • **kwargs (dict) – Other kwargs that will be passed to SALib.

Returns:

si_dct – A dict of Sobol analysis results.

Return type:

dict

Examples

>>> # Import a pre-constructed system model
>>> from qsdsan import stats as s
>>> from qsdsan.utils import create_example_model
>>> model = create_example_model(evaluate=False)
>>> # Sobol analysis requires special samples
>>> inputs = s.define_inputs(model)
>>> # Use a small sample size for demonstrative purpose
>>> samples = s.generate_samples(inputs, kind='Sobol', N=10, calc_second_order=True)
>>> model.load_samples(samples)
>>> model.evaluate()
>>> # Error will be raised if `nan_policy` says so
>>> dct = s.sobol_analysis(
...     model, inputs, seed=554, calc_second_order=True, conf_level=0.95,
...     nan_policy='raise') 
Traceback ...
>>> dct = s.sobol_analysis(
...     model, inputs, seed=554, calc_second_order=True, conf_level=0.95,
...     nan_policy='fill_mean')
>>> # Different types of plots
>>> fig, ax = s.plot_sobol_results(dct, metric=model.metrics[-1], kind='STS1')
>>> fig, ax = s.plot_sobol_results(
...     dct, metric=model.metrics[-1], kind='STS2', plot_in_diagonal='ST')
>>> fig, ax = s.plot_sobol_results(dct, metric=model.metrics[0], kind='all')

See also

qsdsan.stats

SALib API

stats.plot_sobol_results(metric, ax=None, parameters=(), kind='all', annotate_heatmap=False, plot_in_diagonal='', error_bar=True, file='', close_fig=True)

Visualize the results from Sobol analysis as a bar plot and/or heat map. Total (\(S_{Ti}\)) and main (\(S_{1i}\)) effects can be drawn in the bar plot or diagonal of the heat map; second-order interaction effects between two parameters (\(S_{2ij}\)) are shown in the heat map.

Parameters:

  • result_dct (dict) – Result dictionary generated from qsdsan.stats.sobol_analysis()

  • metric (biosteam.Metric) – The metric of interest for the plot.

  • parameters (biosteam.Parameter) – Single or a Iterable of model parameters whose \({\mu^*}\) will be included in the plot. Will be set to all parameters in retult_dct will be used if not provided.

  • kind (str) –

    Which sensitivity index or indices to plot:

    kind

    returned plot

    ST

    total effects (bar)

    S1

    main effects (bar)

    S2

    interaction effects (heat map)

    STS1

    total and main effects (bar)

    STS2

    total and interaction effects (heat map or bar and heat map)

    S1S2

    main and interaction effects (heat map or bar and heat map)

    all

    all effects (bar and heat map)

  • annotate_heatmap (bool) – Whether to annotate the index values in the heat map.

  • plot_in_diagonal (str) – Plot total or main effects in the diagonal of the interaction heat map, can be “ST”, “S1”, or “”. This is applicable when kind is “STS2”, “S1S2”, or “all”.

  • error_bar (bool) – Whether to include the confidence interval as error bars in the plot, this is only applicable for the bar plot.

  • file (str) – If provided, the generated figure will be saved as a png file.

  • close_fig (bool) – Whether to close the figure (if not close, new figure will be overlaid on the current figure).

Returns:

  • figure (matplotlib.figure.Figure) – The generated figure.

  • axis (matplotlib.axes._subplots.AxesSubplot) – The generated figure axis. If generating bar plot and heat map, a tuple of two axes will be returned for the respective plot.

See also

sobol_analysis()

qsdsan.stats

Examples

Bar plot for total and main effects

# Run Sobol analysis, here we use a small sample size for demonstrative purpose
sobol_samples = s.generate_samples(inputs, kind='Sobol', N=10,
                                   calc_second_order=True)

evaluate(modelA, sobol_samples)

dct = s.sobol_analysis(modelA, inputs, metrics=key_metrics, seed=seed,
                       calc_second_order=True, conf_level=0.95,
                       nan_policy='fill_mean')

fig, ax = s.plot_sobol_results(dct, metric=key_metrics[-1], kind='STS1')

fig.subplots_adjust(left=0.4, top=0.95)
../_images/plot_sobol_sts1.png

Heat map for total and second-order effects

fig, ax = s.plot_sobol_results(dct, metric=key_metrics[-1], kind='STS2',
                               plot_in_diagonal='ST')

for label in ax.get_xticklabels():
    label.set_rotation(45)

fig.subplots_adjust(left=0.4, bottom=0.4)
../_images/plot_sobol_sts2.png

Bar plot and heat map for total, main, and second-order effects

fig, ax = s.plot_sobol_results(dct, metric=key_metrics[-1], kind='all')
../_images/plot_sobol_all.png