diff --git a/doc/source/images/kernel_polynomial_method.py.diff b/doc/source/images/kernel_polynomial_method.py.diff
index 9f440dd8bbe363204f3d1b19c00869e9946f0ebf..92c55e20bccdd466c14e37cb939dc1368e3d5fa8 100644
--- a/doc/source/images/kernel_polynomial_method.py.diff
+++ b/doc/source/images/kernel_polynomial_method.py.diff
@@ -56,22 +56,22 @@
 +    )
  
  
- def dos_averaging_example(fsyst):
+ def dos_integrating_example(fsyst):
      spectrum = kwant.kpm.SpectralDensity(fsyst)
  
--    print('identity resolution:', spectrum.average())
+-    print('identity resolution:', spectrum.integrate())
 +    with open('kpm_normalization.txt', 'w') as f:
 +        with redirect_stdout(f):
-+            print('identity resolution:', spectrum.average())
++            print('identity resolution:', spectrum.integrate())
  
      # Fermi energy 0.1 and temperature 0.2
      fermi = lambda E: 1 / (np.exp((E - 0.1) / 0.2) + 1)
      n_states = len(fsyst.sites)  # 1 degree of freedom per site
  
--    print('number of filled states:', n_states * spectrum.average(fermi))
+-    print('number of filled states:', n_states * spectrum.integrate(fermi))
 +    with open('kpm_total_states.txt', 'w') as f:
 +        with redirect_stdout(f):
-+            print('number of filled states:', n_states * spectrum.average(fermi))
++            print('number of filled states:', n_states * spectrum.integrate(fermi))
  
  
  def increasing_accuracy_example(fsyst):
diff --git a/doc/source/tutorial/kernel_polynomial_method.py b/doc/source/tutorial/kernel_polynomial_method.py
index 7e8aa16d2d11f1b3b116ccf76b8d16b0b2a50b34..1f4a52d73a5b44daecd125de263b9140cf8850ee 100644
--- a/doc/source/tutorial/kernel_polynomial_method.py
+++ b/doc/source/tutorial/kernel_polynomial_method.py
@@ -85,20 +85,20 @@ def simple_dos_example():
     ])
 
 
-def dos_averaging_example(fsyst):
+def dos_integrating_example(fsyst):
     spectrum = kwant.kpm.SpectralDensity(fsyst)
 
-#HIDDEN_BEGIN_av1
-    print('identity resolution:', spectrum.average())
-#HIDDEN_END_av1
+#HIDDEN_BEGIN_int1
+    print('identity resolution:', spectrum.integrate())
+#HIDDEN_END_int1
 
-#HIDDEN_BEGIN_av2
+#HIDDEN_BEGIN_int2
     # Fermi energy 0.1 and temperature 0.2
     fermi = lambda E: 1 / (np.exp((E - 0.1) / 0.2) + 1)
     n_states = len(fsyst.sites)  # 1 degree of freedom per site
 
-    print('number of filled states:', n_states * spectrum.average(fermi))
-#HIDDEN_END_av2
+    print('number of filled states:', n_states * spectrum.integrate(fermi))
+#HIDDEN_END_int2
 
 
 def increasing_accuracy_example(fsyst):
@@ -206,7 +206,7 @@ def main():
 
     fsyst = make_syst().finalized()
 
-    dos_averaging_example(fsyst)
+    dos_integrating_example(fsyst)
     increasing_accuracy_example(fsyst)
     operator_example(fsyst)
     ldos_example(fsyst)
diff --git a/doc/source/tutorial/kpm.rst b/doc/source/tutorial/kpm.rst
index f40d1920231ee7602e9d93c2f472a6e9cb023e22..1527c34c28cf345d18f9e0ce426a448878d0225c 100644
--- a/doc/source/tutorial/kpm.rst
+++ b/doc/source/tutorial/kpm.rst
@@ -113,26 +113,26 @@ evaluate the density of states.
 .. image:: ../images/kpm_dos.*
 
 In addition to being called like functions, `~kwant.kpm.SpectralDensity`
-objects also have a method `~kwant.kpm.SpectralDensity.average` which can be
+objects also have a method `~kwant.kpm.SpectralDensity.integrate` which can be
 used to integrate the density of states against some distribution function over
 the whole spectrum. If no distribution function is specified, then the uniform
 distribution is used:
 
 .. literalinclude:: kernel_polynomial_method.py
-    :start-after: #HIDDEN_BEGIN_av1
-    :end-before: #HIDDEN_END_av1
+    :start-after: #HIDDEN_BEGIN_int1
+    :end-before: #HIDDEN_END_int1
 
 .. literalinclude:: ../images/kpm_normalization.txt
 
 We see that the integral of the density of states is normalized to 1. If
-we wish to calculate, say, the average number of states populated in
+we wish to calculate, say, the number of states populated in
 equilibrium, then we should integrate with respect to a Fermi-Dirac
 distribution and multiply by the total number of available states in
 the system:
 
 .. literalinclude:: kernel_polynomial_method.py
-    :start-after: #HIDDEN_BEGIN_av2
-    :end-before: #HIDDEN_END_av2
+    :start-after: #HIDDEN_BEGIN_int2
+    :end-before: #HIDDEN_END_int2
 
 .. literalinclude:: ../images/kpm_total_states.txt
 
diff --git a/kwant/kpm.py b/kwant/kpm.py
index f55b532b9f3cb2923eb94eedd4b5a77cadbf61f2..a4ad191e511a6f7ed7b4d97c83cf7c3157e23436 100644
--- a/kwant/kpm.py
+++ b/kwant/kpm.py
@@ -255,7 +255,7 @@ class SpectralDensity:
             return np.transpose(np.polynomial.chebyshev.chebval(
                     rescaled_energy, coef_cheb) / g_e).real
 
-    def average(self, distribution_function=None):
+    def integrate(self, distribution_function=None):
         """Returns the total spectral density.
 
         Returns the integral over the whole spectrum with an optional
diff --git a/kwant/tests/test_kpm.py b/kwant/tests/test_kpm.py
index 5ea3dbf76d508b9c28f415f223b5ce662cc7d421..3caf1d3d2e7e83b002283a35c898f6e64b27e99f 100644
--- a/kwant/tests/test_kpm.py
+++ b/kwant/tests/test_kpm.py
@@ -290,7 +290,7 @@ def test_kwant_op_current():
     assert_allclose(spectrum_syst.densities, spectrum.densities)
 
 
-def test_kwant_op_average():
+def test_kwant_op_integrate():
     """Check that the kwant.operator.Density gives the same result as the
     identity operator.
     """
@@ -305,8 +305,8 @@ def test_kwant_op_average():
 
     assert spectrum_syst.densities.shape[1] == ham.shape[0]
     # same algorithms are used so these arrays are equal up to TOL
-    assert_allclose(np.sum(spectrum_syst.average(distribution_function=ones)),
-                    spectrum.average())
+    assert_allclose(np.sum(spectrum_syst.integrate(distribution_function=ones)),
+                    spectrum.integrate())
 
 
 # ## test for methods to work as expected
@@ -495,16 +495,16 @@ def test_call():
     # different algorithms are used so these arrays are equal up to TOL_SP
     assert_allclose_sp(densities_array, spectrum.densities)
 
-# ### check average
+# ### check integrate
 
 
-def test_average():
+def test_integrate():
     ham = kwant.rmt.gaussian(dim)
     spectrum = make_spectrum(ham, p)
     ones = lambda x: np.ones_like(x)
-    assert np.abs(spectrum.average() - simps(spectrum.densities,
+    assert np.abs(spectrum.integrate() - simps(spectrum.densities,
                                              x=spectrum.energies)) < TOL_SP
-    assert np.abs(spectrum.average() - spectrum.average(
+    assert np.abs(spectrum.integrate() - spectrum.integrate(
         distribution_function=ones)) < TOL
 
 # ### check increase_energy_resolution