Source code for SALib.test_functions.oakley2004

import numpy as np

from SALib import ProblemSpec


[docs] def evaluate(X: np.ndarray, A: np.ndarray, M: np.ndarray) -> np.ndarray: """Test function taken from Oakley and O'Hagan (2004) (see Eqn. 21 in [1]) References ---------- .. [1] Oakley, J.E., O’Hagan, A., 2004. Probabilistic sensitivity analysis of complex models: a Bayesian approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 66, 751–769. https://doi.org/10.1111/j.1467-9868.2004.05304.x """ a1, a2, a3 = A num_ins = X.shape[0] Y = np.zeros(num_ins) for i in range(num_ins): x_i = X[i] Y[i] = a1 @ x_i + a2 @ np.sin(x_i) + a3 @ np.cos(x_i) + x_i.T @ (M @ x_i) a3.T.dot(np.cos(X[0])) return Y
if __name__ == "__main__": # Raw values taken from: http://www.jeremy-oakley.staff.shef.ac.uk/psa_example.txt M = np.array( [ -2.25e-02, -1.85e-01, 1.34e-01, 3.69e-01, 1.72e-01, 1.37e-01, -4.40e-01, -8.14e-02, 7.13e-01, -4.44e-01, 5.04e-01, -2.41e-02, -4.59e-02, 2.17e-01, 5.59e-02, 2.57e-01, 5.38e-02, 2.58e-01, 2.38e-01, -5.91e-01, -8.16e-02, -2.87e-01, 4.16e-01, 4.98e-01, 8.39e-02, -1.11e-01, 3.32e-02, -1.40e-01, -3.10e-02, -2.23e-01, -5.60e-02, 1.95e-01, 9.55e-02, -2.86e-01, -1.44e-01, 2.24e-01, 1.45e-01, 2.90e-01, 2.31e-01, -3.19e-01, -2.90e-01, -2.10e-01, 4.31e-01, 2.44e-02, 4.49e-02, 6.64e-01, 4.31e-01, 2.99e-01, -1.62e-01, -3.15e-01, -3.90e-01, 1.77e-01, 5.80e-02, 1.72e-01, 1.35e-01, -3.53e-01, 2.51e-01, -1.88e-02, 3.65e-01, -3.25e-01, -1.21e-01, 1.25e-01, 1.07e-01, 4.66e-02, -2.17e-01, 1.95e-01, -6.55e-02, 2.44e-02, -9.68e-02, 1.94e-01, 3.34e-01, 3.13e-01, -8.36e-02, -2.53e-01, 3.73e-01, -2.84e-01, -3.28e-01, -1.05e-01, -2.21e-01, -1.37e-01, -1.44e-01, -1.15e-01, 2.24e-01, -3.04e-02, -5.15e-01, 1.73e-02, 3.90e-02, 3.61e-01, 3.09e-01, 5.00e-02, -7.79e-02, 3.75e-03, 8.87e-01, -2.66e-01, -7.93e-02, -4.27e-02, -1.87e-01, -3.56e-01, -1.75e-01, 8.87e-02, 4.00e-01, -5.60e-02, 1.37e-01, 2.15e-01, -1.13e-02, -9.23e-02, 5.92e-01, 3.13e-02, -3.31e-02, -2.43e-01, -9.98e-02, 3.45e-02, 9.51e-02, -3.38e-01, 6.39e-03, -6.12e-01, 8.13e-02, 8.87e-01, 1.43e-01, 1.48e-01, -1.32e-01, 5.29e-01, 1.27e-01, 4.51e-02, 5.84e-01, 3.73e-01, 1.14e-01, -2.95e-01, -5.70e-01, 4.63e-01, -9.41e-02, 1.40e-01, -3.86e-01, -4.49e-01, -1.46e-01, 5.81e-02, -3.23e-01, 9.31e-02, 7.24e-02, -5.69e-01, 5.26e-01, 2.37e-01, -1.18e-02, 7.18e-02, 7.83e-02, -1.34e-01, 2.27e-01, 1.44e-01, -4.52e-01, -5.56e-01, 6.61e-01, 3.46e-01, 1.41e-01, 5.19e-01, -2.80e-01, -1.60e-01, -6.84e-02, -2.04e-01, 6.97e-02, 2.31e-01, -4.44e-02, -1.65e-01, 2.16e-01, 4.27e-03, -8.74e-02, 3.16e-01, -2.76e-02, 1.34e-01, 1.35e-01, 5.40e-02, -1.74e-01, 1.75e-01, 6.03e-02, -1.79e-01, -3.11e-01, -2.54e-01, 2.58e-02, -4.30e-01, -6.23e-01, -3.40e-02, -2.90e-01, 3.41e-02, 3.49e-02, -1.21e-01, 2.60e-02, -3.35e-01, -4.14e-01, 5.32e-02, -2.71e-01, -2.63e-02, 4.10e-01, 2.66e-01, 1.56e-01, -1.87e-01, 1.99e-02, -2.44e-01, -4.41e-01, 1.26e-02, 2.49e-01, 7.11e-02, 2.46e-01, 1.75e-01, 8.53e-03, 2.51e-01, -1.47e-01, -8.46e-02, 3.69e-01, -3.00e-01, 1.10e-01, -7.57e-01, 4.15e-02, -2.60e-01, 4.64e-01, -3.61e-01, -9.50e-01, -1.65e-01, 3.09e-03, 5.28e-02, 2.25e-01, 3.84e-01, 4.56e-01, -1.86e-01, 8.23e-03, 1.67e-01, 1.60e-01, ] ).reshape(15, 15) A = np.array( [ [ 0.0118, 0.0456, 0.2297, 0.0393, 0.1177, 0.3865, 0.3897, 0.6061, 0.6159, 0.4005, 1.0741, 1.1474, 0.7880, 1.1242, 1.1982, ], [ 0.4341, 0.0887, 0.0512, 0.3233, 0.1489, 1.0360, 0.9892, 0.9672, 0.8977, 0.8083, 1.8426, 2.4712, 2.3946, 2.0045, 2.2621, ], [ 0.1044, 0.2057, 0.0774, 0.2730, 0.1253, 0.7526, 0.8570, 1.0331, 0.8388, 0.7970, 2.2145, 2.0382, 2.4004, 2.0541, 1.9845, ], ] ) sp = ProblemSpec( { "names": ["x{}".format(i) for i in range(1, 16)], "bounds": [ [0.0, 0.835], ] * 15, "dists": ["norm"] * 15, } ) ( sp.sample_latin(2048, seed=101) .evaluate(evaluate, A, M) .analyze_rbd_fast(seed=101, num_resamples=100) ) print(sp) # analytic S1 values analytic = np.array( [ 0.00156, 0.000186, 0.001307, 0.003045, 0.002905, 0.023035, 0.024151, 0.026517, 0.046036, 0.014945, 0.101823, 0.135708, 0.101989, 0.105169, 0.122818, ] ) S1 = sp.analysis.to_df() S1["lower"] = S1["S1"] - S1["S1_conf"] S1["analytic"] = analytic S1["upper"] = S1["S1"] + S1["S1_conf"] print(np.all((analytic >= S1["lower"]) & (analytic <= S1["upper"]))) print((analytic >= S1["lower"]) & (analytic <= S1["upper"]))