<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Conservative Estimation using a Grid Seach Minimization
This notebook illustrates the different steps for a conservative estimation using a grid search minimization.
Classic Libraries
Step1:  Additive model
The first example of conservative estimation consider an additive model $\eta 
Step2:  Dimension 2
We consider the problem in dimension $d=2$ and a number of pairs $p=1$ for gaussian margins.
Step3:  Copula families
We consider a gaussian copula for this first example
Step4:  Estimations
We create an instance of the main class for a conservative estimate.
Step5:  First, we compute the quantile at independence
Step6:  We aim to minimize the output quantile. To do that, we create a q_func object from the function quantile_func to associate a probability $\alpha$ to a function that computes the empirical quantile from a given sample.
Step7:  The computation returns a DependenceResult instance. This object gather the informations of the computation. It also computes the output quantity of interest (which can also be changed).
Step8:  A boostrap can be done on the output quantity
Step9:  And we can plot it
Step10:  Grid Search Approach
Firstly, we consider a grid search approach in order to compare the perfomance with the iterative algorithm. The discretization can be made on the parameter space or on other concordance measure such as the kendall's Tau. This below example shows a grid-search on the parameter space.
Step11:  The computation returns a ListDependenceResult which is a list of DependenceResult instances and some bonuses.
Step12:  Lets set the quantity function and search for the minimum among the grid results.
Step13:  We can plot the result in grid results. The below figure shows the output quantiles in function of the dependence parameters.
Step14:  As for the individual problem, we can do a boostrap also, for each parameters. Because we have $K$ parameters, we can do a bootstrap for the $K$ samples, compute the $K$ quantiles for all the bootstrap and get the minimum quantile for each bootstrap.
Step15:  For the parameter that have the most occurence for the minimum, we compute its bootstrap mean.
Step16:  Kendall's Tau
An interesting feature is to convert the dependence parameters to Kendall's Tau values.
Step17:  As we can see, the bounds
With bounds on the dependencies
An interesting option in the ConservativeEstimate class is to bound the dependencies, due to some prior informations.
Step18:  Saving the results
It is usefull to save the result in a file to load it later and compute other quantities or anything you need!
Step19:  Taking the extreme values of the dependence parameter
If the output quantity of interest seems to have a monotonicity with the dependence parameter, it is better to directly take the bounds of the dependence problem. Obviously, the minimum should be at the edges of the design space
Step20:  Higher Dimension
We consider the problem in dimension $d=5$.
Step21:  Copula families with one dependent pair
We consider a gaussian copula for this first example, but for the moment only one pair is dependent.
Step22:  We reset the families and bounds for the current instance. (I don't want to create a new instance, just to check if the setters are good).
Step23:  Let's do the grid search to see
Step24:  The quantile is lower compare to the problem of dimension 1. Indeed, there is more variables, more uncertainty, so a larger deviation of the output.
Step25:  Copula families with all dependent pairs
We consider a gaussian copula for this first example, but for the moment only one pair is dependent.
Step26:  With one fixed pair
Step27:  Save the used grid and load it again
Step28:  Then gather the results from the same grid with the same configurations
Step29:  Because the configurations are the same, we can gather the results from two different runs
<ASSISTANT_TASK:>
Python Code:
import openturns as ot
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline
%load_ext autoreload
%autoreload 2
random_state = 123
np.random.seed(random_state)
Explanation: Conservative Estimation using a Grid Seach Minimization
This notebook illustrates the different steps for a conservative estimation using a grid search minimization.
Classic Libraries
End of explanation
from depimpact.tests import func_sum
help(func_sum)
Explanation: Additive model
The first example of conservative estimation consider an additive model $\eta : \mathbb R^d \rightarrow \mathbb R$ with Gaussian margins. The objectives are to estimate a quantity of interest $\mathcal C(Y)$ of the model output distribution. Unfortunately, the dependence structure is unknown. In order to be conservative we aim to give bounds to $\mathcal C(Y)$.
The model
This example consider the simple additive example.
End of explanation
dim = 2
margins = [ot.Normal()]*dim
Explanation: Dimension 2
We consider the problem in dimension $d=2$ and a number of pairs $p=1$ for gaussian margins.
End of explanation
families = np.zeros((dim, dim), dtype=int)
families[1, 0] = 1
Explanation: Copula families
We consider a gaussian copula for this first example
End of explanation
from depimpact import ConservativeEstimate
quant_estimate = ConservativeEstimate(model_func=func_sum, margins=margins, families=families)
Explanation: Estimations
We create an instance of the main class for a conservative estimate.
End of explanation
n = 1000
indep_result = quant_estimate.independence(n_input_sample=n, random_state=random_state)
Explanation: First, we compute the quantile at independence
End of explanation
from depimpact import quantile_func
alpha = 0.05
q_func = quantile_func(alpha)
indep_result.q_func = q_func
Explanation: We aim to minimize the output quantile. To do that, we create a q_func object from the function quantile_func to associate a probability $\alpha$ to a function that computes the empirical quantile from a given sample.
End of explanation
sns.jointplot(indep_result.input_sample[:, 0], indep_result.input_sample[:, 1]);
h = sns.distplot(indep_result.output_sample_id, axlabel='Model output', label="Output Distribution")
plt.plot([indep_result.quantity]*2, h.get_ylim(), label='Quantile at %d%%' % (alpha*100))
plt.legend(loc=0)
print('Output quantile :', indep_result.quantity)
Explanation: The computation returns a DependenceResult instance. This object gather the informations of the computation. It also computes the output quantity of interest (which can also be changed).
End of explanation
indep_result.compute_bootstrap(n_bootstrap=5000)
Explanation: A boostrap can be done on the output quantity
End of explanation
sns.distplot(indep_result.bootstrap_sample, axlabel='Output quantile');
ci = [0.025, 0.975]
quantity_ci = indep_result.compute_quantity_bootstrap_ci(ci)
h = sns.distplot(indep_result.output_sample_id, axlabel='Model output', label="Output Distribution")
plt.plot([indep_result.quantity]*2, h.get_ylim(), 'g-', label='Quantile at %d%%' % (alpha*100))
plt.plot([quantity_ci[0]]*2, h.get_ylim(), 'g--', label='%d%% confidence intervals' % ((1. - (ci[0] + 1. - ci[1]))*100))
plt.plot([quantity_ci[1]]*2, h.get_ylim(), 'g--')
plt.legend(loc=0)
print('Quantile at independence: %.2f with a C.O.V at %.1f %%' % (indep_result.boot_mean, indep_result.boot_cov))
Explanation: And we can plot it
End of explanation
K = 20
n = 10000
grid_type = 'lhs'
dep_measure = 'parameter'
grid_result = quant_estimate.gridsearch(n_dep_param=K, n_input_sample=n, grid_type=grid_type, dep_measure=dep_measure, 
                                        random_state=random_state)
Explanation: Grid Search Approach
Firstly, we consider a grid search approach in order to compare the perfomance with the iterative algorithm. The discretization can be made on the parameter space or on other concordance measure such as the kendall's Tau. This below example shows a grid-search on the parameter space.
End of explanation
print('The computation did %d model evaluations.' % (grid_result.n_evals))
Explanation: The computation returns a ListDependenceResult which is a list of DependenceResult instances and some bonuses.
End of explanation
grid_result.q_func = q_func
min_result = grid_result.min_result
print('Minimum quantile: {} at param: {}'.format(min_result.quantity, min_result.dep_param))
Explanation: Lets set the quantity function and search for the minimum among the grid results.
End of explanation
plt.plot(grid_result.dep_params, grid_result.quantities, '.', label='Quantiles')
plt.plot(min_result.dep_param[0], min_result.quantity, 'ro', label='minimum')
plt.xlabel('Dependence parameter')
plt.ylabel('Quantile value')
plt.legend(loc=0);
Explanation: We can plot the result in grid results. The below figure shows the output quantiles in function of the dependence parameters.
End of explanation
grid_result.compute_bootstraps(n_bootstrap=500)
boot_min_quantiles = grid_result.bootstrap_samples.min(axis=0)
boot_argmin_quantiles = grid_result.bootstrap_samples.argmin(axis=0).ravel().tolist()
boot_min_params = [grid_result.dep_params[idx][0] for idx in boot_argmin_quantiles]
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
sns.distplot(boot_min_quantiles, axlabel="Minimum quantiles", ax=axes[0])
sns.distplot(boot_min_params, axlabel="Parameters of the minimum", ax=axes[1])
Explanation: As for the individual problem, we can do a boostrap also, for each parameters. Because we have $K$ parameters, we can do a bootstrap for the $K$ samples, compute the $K$ quantiles for all the bootstrap and get the minimum quantile for each bootstrap.
End of explanation
 # The parameter with most occurence
boot_id_min = max(set(boot_argmin_quantiles), key=boot_argmin_quantiles.count)
boot_min_result = grid_result[boot_id_min]
boot_mean = boot_min_result.bootstrap_sample.mean()
boot_std = boot_min_result.bootstrap_sample.std()
print('Worst Quantile: {} at {} with a C.O.V of {} %'.format(boot_min_result.boot_mean, min_result.dep_param, boot_min_result.boot_cov*100.))
Explanation: For the parameter that have the most occurence for the minimum, we compute its bootstrap mean.
End of explanation
plt.plot(grid_result.kendalls, grid_result.quantities, '.', label='Quantiles')
plt.plot(min_result.kendall_tau, min_result.quantity, 'ro', label='Minimum quantile')
plt.xlabel("Kendall's tau")
plt.ylabel('Quantile')
plt.legend(loc=0);
Explanation: Kendall's Tau
An interesting feature is to convert the dependence parameters to Kendall's Tau values.
End of explanation
bounds_tau = np.asarray([[0., 0.7], [0.1, 0.]])
quant_estimate.bounds_tau = bounds_tau
K = 20
n = 10000
grid_type = 'lhs'
grid_result = quant_estimate.gridsearch(n_dep_param=K, n_input_sample=n, grid_type=grid_type, random_state=random_state)
grid_result.q_func = q_func
min_result = grid_result.min_result
print('Minimum quantile: {} at param: {}'.format(min_result.quantity, min_result.dep_param))
plt.plot(grid_result.dep_params, grid_result.quantities, '.', label='Quantiles')
plt.plot(min_result.dep_param[0], min_result.quantity, 'ro', label='minimum')
plt.xlabel('Dependence parameter')
plt.ylabel('Quantile value')
plt.legend(loc=0);
Explanation: As we can see, the bounds
With bounds on the dependencies
An interesting option in the ConservativeEstimate class is to bound the dependencies, due to some prior informations.
End of explanation
filename = './result.hdf'
grid_result.to_hdf(filename)
from dependence import ListDependenceResult
load_grid_result = ListDependenceResult.from_hdf(filename, q_func=q_func, with_input_sample=False)
np.testing.assert_array_equal(grid_result.output_samples, load_grid_result.output_samples)
import os
os.remove(filename)
Explanation: Saving the results
It is usefull to save the result in a file to load it later and compute other quantities or anything you need!
End of explanation
K = None
n = 1000
grid_type = 'vertices'
grid_result = quant_estimate.gridsearch(n_dep_param=K, n_input_sample=n, grid_type=grid_type, random_state=random_state)
grid_result.q_func = q_func
print("Kendall's Tau : {}, Quantile: {}".format(grid_result.kendalls.ravel(), grid_result.quantities))
from depimpact.plots import matrix_plot_input
matrix_plot_input(grid_result.min_result);
Explanation: Taking the extreme values of the dependence parameter
If the output quantity of interest seems to have a monotonicity with the dependence parameter, it is better to directly take the bounds of the dependence problem. Obviously, the minimum should be at the edges of the design space
End of explanation
dim = 5
quant_estimate.margins = [ot.Normal()]*dim
Explanation: Higher Dimension
We consider the problem in dimension $d=5$.
End of explanation
families = np.zeros((dim, dim), dtype=int)
families[2, 0] = 1
quant_estimate.families = families
families
quant_estimate.bounds_tau = None
quant_estimate.bounds_tau
Explanation: Copula families with one dependent pair
We consider a gaussian copula for this first example, but for the moment only one pair is dependent.
End of explanation
quant_estimate.vine_structure
Explanation: We reset the families and bounds for the current instance. (I don't want to create a new instance, just to check if the setters are good).
End of explanation
K = 20
n = 10000
grid_type = 'vertices'
grid_result = quant_estimate.gridsearch(n_dep_param=K, n_input_sample=n, grid_type=grid_type, random_state=random_state)
Explanation: Let's do the grid search to see
End of explanation
grid_result.q_func = q_func
min_result = grid_result.min_result
print('Worst Quantile: {} at {}'.format(min_result.quantity, min_result.dep_param))
matrix_plot_input(min_result)
plt.plot(grid_result.dep_params, grid_result.quantities, '.', label='Quantiles')
plt.plot(min_result.dep_param[0], min_result.quantity, 'ro', label='Minimum')
plt.xlabel('Dependence parameter')
plt.ylabel('Quantile value')
plt.legend(loc=0);
Explanation: The quantile is lower compare to the problem of dimension 1. Indeed, there is more variables, more uncertainty, so a larger deviation of the output.
End of explanation
families = np.zeros((dim, dim), dtype=int)
for i in range(1, dim):
    for j in range(i):
        families[i, j] = 1
quant_estimate.margins = margins
quant_estimate.families = families
quant_estimate.vine_structure = None
quant_estimate.bounds_tau = None
quant_estimate.bounds_tau
K = 100
n = 1000
grid_type = 'lhs'
grid_result = quant_estimate.gridsearch(n_dep_param=K, n_input_sample=n, grid_type=grid_type, random_state=random_state)
min_result = grid_result.min_result
print('Worst Quantile: {0} at {1}'.format(min_result.quantity, min_result.dep_param))
Explanation: Copula families with all dependent pairs
We consider a gaussian copula for this first example, but for the moment only one pair is dependent.
End of explanation
families[3, 2] = 0
quant_estimate = ConservativeEstimate(model_func=func_sum, margins=margins, families=families)
K = 100
n = 10000
grid_type = 'lhs'
grid_result = quant_estimate.gridsearch(n_dep_param=K, n_input_sample=n, grid_type=grid_type, 
                                                 q_func=q_func, random_state=random_state)
min_result = grid_result.min_result
print('Worst Quantile: {0} at {1}'.format(min_result.quantity, min_result.dep_param))
grid_result.vine_structure
from depimpact.plots import matrix_plot_input
matrix_plot_input(min_result)
Explanation: With one fixed pair
End of explanation
K = 100
n = 1000
grid_type = 'lhs'
grid_result_1 = quant_estimate.gridsearch(n_dep_param=K, n_input_sample=n, grid_type=grid_type, save_grid=True, grid_path='./output')
grid_result_2 = quant_estimate.gridsearch(n_dep_param=K, n_input_sample=n, grid_type=grid_type, 
                                                 q_func=q_func, use_grid=0, grid_path='./output')
Explanation: Save the used grid and load it again
End of explanation
grid_result_1.n_input_sample, grid_result_2.n_input_sample
grid_result = grid_result_1 + grid_result_2
Explanation: Then gather the results from the same grid with the same configurations
End of explanation
grid_result.n_input_sample
Explanation: Because the configurations are the same, we can gather the results from two different runs
End of explanation
<END_TASK>