Merge pull request #46 from willu47/feat_mu_conf_speedup

Vectorised mu_star_confidence

SALib/analyze/morris.py

@@ -140,23 +140,25 @@ def compute_elementary_effects(model_inputs, model_outputs, trajectory_size, del
 
 
 def compute_mu_star_confidence(ee, num_trajectories, num_resamples, conf_level):
-
+    '''
+    Uses bootstrapping where the elementary effects are resampled with replacement
+    to produce a histogram of resampled mu_star metrics.
+    This resample is used to produce a confidence interval.
+    '''
     ee_resampled = np.empty([num_trajectories])
     mu_star_resampled = np.empty([num_resamples])
 
     if conf_level < 0 or conf_level > 1:
         raise ValueError("Confidence level must be between 0-1.")
 
-    for i in range(num_resamples):
-        for j in range(num_trajectories):
-
-            index = np.random.randint(0, num_trajectories)
-            ee_resampled[j] = ee[index]
-
-        mu_star_resampled[i] = np.average(np.abs(ee_resampled))
+    resample_index = np.floor(np.random.rand(num_resamples * num_trajectories)*len(ee)).astype(int).reshape(num_resamples, num_trajectories)
+    ee_resampled = ee[resample_index]
+    # Compute average of the absolute values over each of the resamples
+    mu_star_resampled = np.average(np.abs(ee_resampled), axis=1)
 
     return norm.ppf(0.5 + conf_level / 2) * mu_star_resampled.std(ddof=1)
 
+
 if __name__ == "__main__":
 
     parser = common_args.create()

SALib/tests/test_analyze_morris.py

@@ -12,8 +12,25 @@
     compute_grouped_mu_star
 
 
-def test_analysis_of_morris_results():
+def test_compute_mu_star_confidence():
+    '''
+    Tests that compute mu_star_confidence is computed correctly
+    '''
+
+    ee = np.array([2.52, 2.01, 2.30, 0.66, 0.93, 1.3], dtype=np.float)
+    num_trajectories = 6
+    num_resamples = 1000
+    conf_level = 0.95
+
+    actual = compute_mu_star_confidence(ee, num_trajectories, num_resamples, conf_level)
+    expected = 0.5
+    assert_allclose(actual, expected, atol=1e-01)
 
+
+def test_analysis_of_morris_results():
+    '''
+    Tests a one-dimensional vector of results
+    '''
     model_input = np.array([[0, 1. / 3], [0, 1],       [2. / 3, 1],
                              [0, 1. / 3],   [2. / 3, 1. / 3], [2. / 3, 1],
                              [2. / 3, 0],   [2. / 3, 2. / 3], [0, 2. / 3],
@@ -60,6 +77,12 @@ def test_conf_level_within_zero_one_bounds():
 
 
 def test_compute_elementary_effects():
+    '''
+    Inputs for elementary effects taken from Exercise 5 from Saltelli (2008).
+    See page 140-145.
+    `model_inputs` are from trajectory t_1 from table 3.10 on page 141.
+    `desired` is equivalent to column t_1 in table 3.12 on page 145.
+    '''
     model_inputs = np.array([
                             [1.64, -1.64, -1.64, 0.39, -0.39, 0.39, -1.64, -
                                 1.64, -0.39, -0.39, 1.64, 1.64, -0.39, 0.39, 1.64],