diff --git a/python/damask/util.py b/python/damask/util.py
index 43e86314e..12d4fd6cc 100644
--- a/python/damask/util.py
+++ b/python/damask/util.py
@@ -20,6 +20,7 @@ __all__=[
          'execute',
          'show_progress',
          'scale_to_coprime',
+         'hybrid_IA',
          'return_message',
          'extendableOption',
          'execution_stamp'
@@ -187,6 +188,23 @@ def execution_stamp(class_name,function_name=None):
     return f'damask.{class_name}{_function_name} v{version} ({now})'
 
 
+def hybrid_IA(dist,N,seed=None):
+    rng = np.random.default_rng(seed)
+    N_opt_samples = max(np.count_nonzero(dist),N)                                                   # random subsampling if too little samples requested
+
+    scale_,scale,inc_factor = (0.0,float(N_opt_samples),1.0)
+    repeats = np.rint(scale*dist).astype(int)
+    N_inv_samples = np.sum(repeats)
+    while (not np.isclose(scale, scale_)) and (N_inv_samples != N_opt_samples):
+        scale_,scale,inc_factor = (scale,scale+inc_factor*0.5*(scale - scale_), inc_factor*2.0) \
+                                   if N_inv_samples < N_opt_samples else \
+                                  (scale_,0.5*(scale_ + scale), 1.0)
+        repeats = np.rint(scale*dist).astype(int)
+        N_inv_samples = np.sum(repeats)
+
+    return np.repeat(np.arange(len(dist)),repeats)[rng.permutation(N_inv_samples)[:N]]
+
+
 ####################################################################################################
 # Classes
 ####################################################################################################
diff --git a/python/tests/test_util.py b/python/tests/test_util.py
index 053045741..0b5593e8e 100644
--- a/python/tests/test_util.py
+++ b/python/tests/test_util.py
@@ -1,5 +1,7 @@
 import pytest
 import numpy as np
+from scipy import stats
+
 from damask import util
 
 
@@ -31,3 +33,14 @@ class TestUtil:
     def test_lackofprecision(self):
         with pytest.raises(ValueError):
             util.scale_to_coprime(np.array([1/333.333,1,1]))
+
+
+    @pytest.mark.parametrize('rv',[stats.rayleigh(),stats.weibull_min(1.2),stats.halfnorm(),stats.pareto(2.62)])
+    def test_hybridIA(self,rv):
+        bins = np.linspace(0,10,100000)
+        centers = (bins[1:]+bins[:-1])/2
+        N_samples = bins.shape[0]-1000
+        dist = rv.pdf(centers)
+        selected = util.hybrid_IA(dist,N_samples)
+        dist_sampled = np.histogram(centers[selected],bins)[0]/N_samples*np.sum(dist)
+        assert np.sqrt(((dist - dist_sampled) ** 2).mean()) < .025 and selected.shape[0]==N_samples