In [1]:
import openturns as ot
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
%load_ext autoreload
%autoreload 2
random_state = 123
np.random.seed(random_state)
In [2]:
from depimpact.tests import func_overflow, margins_overflow, var_names_overflow, func_sum
from depimpact.tests.test_functions import func_prod
from dependence import iterative_vine_minimize
test_func = func_sum
In [3]:
if test_func == func_overflow:
margins = margins_overflow
dim = len(margins)
else:
dim = 8
margins = [ot.Normal()]*dim
We chose the coefficients of the variables throught the additive function.
In [4]:
if test_func in [func_sum, func_prod]:
coeficients = np.logspace(0., 3., dim+1, endpoint=False)[1:]
coeficients /= coeficients.max()
#coeficients = np.asarray(range(dim))
new_test_func = lambda x: test_func(x, coeficients)
n_plot = 10000
x = np.asarray(ot.ComposedDistribution(margins).getSample(n_plot))
y = test_func(x, coeficients)
fig, axes = plt.subplots(dim, 1, sharex=False, sharey=True, figsize=(4, 2*dim))
for i in range(dim):
ax = axes[i]
ax.plot(x[:, i], y, '.')
ax.set_xlabel(r'$X_{%d}$' % (i+1), fontsize=12)
ax.set_ylabel(r'$y$', fontsize=12)
fig.tight_layout()
else:
new_test_func = test_func
In [5]:
families = np.zeros((dim, dim), dtype=int)
for i in range(1, dim):
for j in range(i):
families[i, j] = 1
In [6]:
from depimpact import ConservativeEstimate, quantile_func
alpha = 0.95
if alpha > 0.5: # Maximizing the quantile
def q_func(x, axis=1):
return - quantile_func(alpha)(x, axis=axis)
else: # Minimizing
q_func = quantile_func(alpha)
quant_estimate = ConservativeEstimate(model_func=new_test_func, margins=margins, families=families)
First, we compute the quantile at independence
In [7]:
n = 10000
indep_result = quant_estimate.independence(n_input_sample=n, q_func=q_func, random_state=random_state)
In [8]:
indep_result.compute_bootstrap(1000)
boot_std = indep_result.bootstrap_sample.std()
boot_mean = indep_result.bootstrap_sample.mean()
print('Quantile at independence: %.2f with a C.O.V at %.1f %%' % (boot_mean, abs(boot_std/boot_mean)*100.))
In [215]:
algorithm_parameters = {
"n_input_sample": 10000,
"n_dep_param_init": 10,
"max_n_pairs": 1,
"grid_type": 'lhs',
"q_func": q_func,
"n_add_pairs": 1,
"n_remove_pairs": 0,
"adapt_vine_structure": True,
"with_bootstrap": False,
"verbose": False,
"iterative_save": False,
"iterative_load": False,
"load_input_samples": False,
"keep_input_samples": False
}
quant_estimate = ConservativeEstimate(model_func=new_test_func, margins=margins, families=families)
iterative_results = iterative_vine_minimize(estimate_object=quant_estimate, **algorithm_parameters)
In [216]:
from depimpact.dependence_plot import matrix_plot_quantities
results = iterative_results[0]
matrix_plot_quantities(results, figsize=(18, 15))
# plt.savefig('output/matrix_plot.png')
In [20]:
algorithm_parameters = {
"n_input_sample": 10000,
"n_dep_param_init": 20,
"max_n_pairs": 3,
"grid_type": 'vertices',
"q_func": q_func,
"n_add_pairs": 2,
"n_remove_pairs": 7,
"adapt_vine_structure": False,
"with_bootstrap": False,
"verbose": True,
"iterative_save": False,
"iterative_load": False,
"load_input_samples": False,
"keep_input_samples": False,
"input_names": var_names_overflow
}
quant_estimate = ConservativeEstimate(model_func=new_test_func, margins=margins, families=families)
iterative_results = iterative_vine_minimize(estimate_object=quant_estimate, **algorithm_parameters)
In [17]:
from depimpact.dependence_plot import plot_iterative_results, matrix_plot_input
plot_iterative_results(iterative_results, indep_result=indep_result, q_func=q_func)
In [18]:
n = 10000
K = 300
quant_estimate.vine_structure = None
grid_results = quant_estimate.gridsearch(n_dep_param=K, n_input_sample=n, q_func=q_func, grid_type='vertices')
In [19]:
min_result = grid_results.min_result
print(min_result.quantity)
In [ ]:
from depimpact.dependence_plot import matrix_plot_input
matrix_plot_input(min_result, margins=margins)
In [46]:
from depimpact.utils import get_pair_id, get_pairs_by_levels, get_possible_structures
pairs_iter = [[1, 0]]
pairs_iter_id = [get_pair_id(dim, pair, with_plus=False) for pair in pairs_iter]
pairs_by_levels = get_pairs_by_levels(dim, pairs_iter_id)
quant_estimate.vine_structure = get_possible_structures(dim, pairs_by_levels)[0]
quant_estimate.vine_structure
Out[46]:
In [47]:
grid_results = quant_estimate.gridsearch(n_dep_param=K, n_input_sample=n, q_func=q_func)
In [48]:
min_result = grid_results.min_result
print(min_result.quantity)
In [49]:
from depimpact.dependence_plot import matrix_plot_input
matrix_plot_input(min_result, margins=margins)
Out[49]:
In [ ]: