diff --git a/exspy/components/eels_arctan.py b/exspy/components/eels_arctan.py
index 99b190f5a..5a52432bf 100644
--- a/exspy/components/eels_arctan.py
+++ b/exspy/components/eels_arctan.py
@@ -17,10 +17,10 @@
 # along with eXSpy. If not, see <https://www.gnu.org/licenses/#GPL>.
 
 
-from hyperspy._components.expression import Expression
+import hyperspy.api as hs
 
 
-class EELSArctan(Expression):
+class EELSArctan(hs.model.components1D.Expression):
     r"""Arctan function component for EELS (with minimum at zero).
 
     .. math::
diff --git a/exspy/components/eels_double_power_law.py b/exspy/components/eels_double_power_law.py
index dd8cb5336..7db5f49c2 100644
--- a/exspy/components/eels_double_power_law.py
+++ b/exspy/components/eels_double_power_law.py
@@ -20,10 +20,10 @@
 import numpy as np
 
 from hyperspy.docstrings.parameters import FUNCTION_ND_DOCSTRING
-from hyperspy._components.expression import Expression
+import hyperspy.api as hs
 
 
-class DoublePowerLaw(Expression):
+class DoublePowerLaw(hs.model.components1D.Expression):
     r"""Double power law component for EELS spectra.
 
     .. math::
diff --git a/exspy/components/eels_vignetting.py b/exspy/components/eels_vignetting.py
index f52266fe2..1b887fc2d 100644
--- a/exspy/components/eels_vignetting.py
+++ b/exspy/components/eels_vignetting.py
@@ -19,7 +19,7 @@
 
 import numpy as np
 from hyperspy.component import Component
-from hyperspy._components.gaussian import Gaussian
+import hyperspy.api as hs
 
 
 class Vignetting(Component):
@@ -46,7 +46,7 @@ def __init__(self):
         self.right.value = np.nan
         self.side_vignetting = False
         self.fix_side_vignetting()
-        self.gaussian = Gaussian()
+        self.gaussian = hs.model.components1D.Gaussian()
         self.gaussian.centre.free, self.gaussian.A.free = False, False
         self.sigma.value = 1.0
         self.gaussian.A.value = 1.0
diff --git a/exspy/components/pes_see.py b/exspy/components/pes_see.py
index 68445ad55..b7bf6b2a6 100644
--- a/exspy/components/pes_see.py
+++ b/exspy/components/pes_see.py
@@ -20,12 +20,13 @@
 import numpy as np
 import logging
 
-from hyperspy._components.expression import Expression
+import hyperspy.api as hs
+
 
 _logger = logging.getLogger(__name__)
 
 
-class SEE(Expression):
+class SEE(hs.model.components1D.Expression):
     r"""Secondary electron emission component for Photoemission Spectroscopy.
 
     .. math::
diff --git a/exspy/misc/eels/base_gos.py b/exspy/misc/eels/base_gos.py
index 49f045d00..58c0a6768 100644
--- a/exspy/misc/eels/base_gos.py
+++ b/exspy/misc/eels/base_gos.py
@@ -144,7 +144,7 @@ def integrateq(self, onset_energy, angle, E0):
             # Perform the integration in a log grid
             qaxis, gos = self.get_qaxis_and_gos(i, qmin, qmax)
             logsqa0qaxis = np.log((a0 * qaxis) ** 2)
-            qint[i] = integrate.simps(gos, logsqa0qaxis)
+            qint[i] = integrate.simpson(gos, x=logsqa0qaxis)
         E = self.energy_axis + energy_shift
         # Energy differential cross section in (barn/eV/atom)
         qint *= (4.0 * np.pi * a0**2.0 * R**2 / E / T * self.subshell_factor) * 1e28
diff --git a/exspy/misc/eels/tools.py b/exspy/misc/eels/tools.py
index 276c87f9c..e7e255467 100644
--- a/exspy/misc/eels/tools.py
+++ b/exspy/misc/eels/tools.py
@@ -60,16 +60,14 @@ def _estimate_gain(
         variance2fit = variance
         average2fit = average
 
-    fit = np.polyfit(average2fit, variance2fit, pol_order)
+    fit = np.polynomial.Polynomial.fit(average2fit, variance2fit, pol_order)
     if weighted is True:
-        from hyperspy._signals.signal1D import Signal1D
-        from hyperspy.models.model1d import Model1D
-        from hyperspy.components1d import Line
+        import hyperspy.api as hs
 
-        s = Signal1D(variance2fit)
+        s = hs.signals.Signal1D(variance2fit)
         s.axes_manager.signal_axes[0].axis = average2fit
-        m = Model1D(s)
-        line = Line()
+        m = s.create_model()
+        line = hs.model.components1D.Polynomial()
         line.a.value = fit[1]
         line.b.value = fit[0]
         m.append(line)
@@ -109,29 +107,30 @@ def estimate_variance_parameters(
     """Find the scale and offset of the Poissonian noise
 
     By comparing an SI with its denoised version (i.e. by PCA),
-    this plots an
-    estimation of the variance as a function of the number of counts
-    and fits a
-    polynomy to the result.
+    this plots an estimation of the variance as a function of the number of counts
+    and fits a polynomial to the result.
 
     Parameters
     ----------
-    noisy_SI, clean_SI : signal1D.Signal1D instances
-    mask : numpy bool array
+    noisy_SI, clean_SI : hyperspy.api.signals.Signal1D
+    mask : numpy.ndarray
         To define the channels that will be used in the calculation.
     pol_order : int
-        The order of the polynomy.
+        The order of the polynomial.
     higher_than: float
         To restrict the fit to counts over the given value.
-    return_results : Bool
-    plot_results : Bool
+    return_results : bool
+        Whether to return the results or not.
+    plot_results : bool
+        Whether to plot the results or not.
     store_results: {True, False, "ask"}, default "ask"
         If True, it stores the result in the signal metadata
 
     Returns
     -------
-    Dictionary with the result of a linear fit to estimate the offset
-    and scale factor
+    dict
+        Dictionary with the result of a linear fit to estimate the offset
+        and scale factor
 
     """
     with noisy_signal.unfolded(), clean_signal.unfolded():
diff --git a/exspy/signals/dielectric_function.py b/exspy/signals/dielectric_function.py
index 047ea90e1..0997aeb68 100644
--- a/exspy/signals/dielectric_function.py
+++ b/exspy/signals/dielectric_function.py
@@ -18,7 +18,7 @@
 
 import numpy as np
 from scipy import constants
-from scipy.integrate import simps, cumtrapz
+from scipy.integrate import simpson, cumtrapz
 
 from hyperspy._signals.complex_signal1d import (
     ComplexSignal1D,
@@ -82,12 +82,12 @@ def get_number_of_effective_electrons(self, nat, cumulative=False):
 
         axis = self.axes_manager.signal_axes[0]
         if cumulative is False:
-            dneff1 = k * simps(
+            dneff1 = k * simpson(
                 (-1.0 / self.data).imag * axis.axis,
                 x=axis.axis,
                 axis=axis.index_in_array,
             )
-            dneff2 = k * simps(
+            dneff2 = k * simpson(
                 self.data.imag * axis.axis, x=axis.axis, axis=axis.index_in_array
             )
             neff1 = self._get_navigation_signal(data=dneff1)
diff --git a/exspy/signals/eds.py b/exspy/signals/eds.py
index 2486cc07c..6a7115870 100644
--- a/exspy/signals/eds.py
+++ b/exspy/signals/eds.py
@@ -24,16 +24,17 @@
 from collections.abc import Iterable
 from matplotlib import pyplot as plt
 
+import hyperspy.api as hs
 from hyperspy import utils
 from hyperspy.signal import BaseSignal
 from hyperspy._signals.signal1d import Signal1D, LazySignal1D
 from exspy.misc.elements import elements as elements_db
 from exspy.misc.eds import utils as utils_eds
 from hyperspy.misc.utils import isiterable
-from hyperspy.utils.markers import Texts, VerticalLines, Lines
 from hyperspy.docstrings.plot import BASE_PLOT_DOCSTRING_PARAMETERS, PLOT1D_DOCSTRING
 from hyperspy.docstrings.signal import LAZYSIGNAL_DOC
 
+
 _logger = logging.getLogger(__name__)
 
 
@@ -1012,14 +1013,14 @@ def _add_vertical_lines_groups(self, position, render_figure=True, **kwargs):
             The position on the signal axis. Each row corresponds to a
             group.
         kwargs
-            keywords argument for :py:class:`~.api.plot.markers.VerticalLine`
+            keywords argument for :class:`hyperspy.api.plot.markers.VerticalLine`
         """
         colors = itertools.cycle(
             np.sort(plt.rcParams["axes.prop_cycle"].by_key()["color"])
         )
 
         for x, color in zip(position, colors):
-            line = VerticalLines(offsets=x, color=color, **kwargs)
+            line = hs.plot.markers.VerticalLines(offsets=x, color=color, **kwargs)
             self.add_marker(line, render_figure=False)
         if render_figure:
             self._render_figure(plot=["signal_plot"])
@@ -1050,12 +1051,12 @@ def add_xray_lines_markers(self, xray_lines, render_figure=True):
                 % utils_eds._get_element_and_line(xray_line)
             )
 
-        line_markers = Lines(
+        line_markers = hs.plot.markers.Lines(
             segments=segments,
             transform="relative",
             color="black",
         )
-        text_markers = Texts(
+        text_markers = hs.plot.markers.Texts(
             offsets=offsets,
             texts=line_names,
             offset_transform="relative",
@@ -1137,7 +1138,7 @@ def _add_background_windows_markers(self, windows_position, render_figure=True):
             x2 = (bw[2] + bw[3]) / 2.0
             segments.append([[x1, y1[0]], [x2, y2[0]]])
         segments = np.array(segments)
-        lines = Lines(segments=segments, color="black")
+        lines = hs.plot.markers.Lines(segments=segments, color="black")
         self.add_marker(lines, render_figure=False)
         if render_figure:
             self._render_figure(plot=["signal_plot"])
diff --git a/exspy/signals/eels.py b/exspy/signals/eels.py
index cb5286b6a..b711df9cf 100644
--- a/exspy/signals/eels.py
+++ b/exspy/signals/eels.py
@@ -28,13 +28,10 @@
 import hyperspy.api as hs
 from hyperspy.signal import BaseSetMetadataItems, BaseSignal
 from hyperspy._signals.signal1d import Signal1D, LazySignal1D
-import hyperspy.axes
-from hyperspy.components1d import PowerLaw
 from hyperspy.misc.utils import display, isiterable, underline
 from hyperspy.misc.math_tools import optimal_fft_size
 
 from hyperspy.ui_registry import add_gui_method, DISPLAY_DT, TOOLKIT_DT
-from hyperspy.utils.markers import Texts, Lines
 from hyperspy.docstrings.signal1d import (
     CROP_PARAMETER_DOC,
     SPIKES_DIAGNOSIS_DOCSTRING,
@@ -622,10 +619,10 @@ def estimating_function(data, threshold=None):
                     if binned:
                         return data.sum()
                     else:
-                        from scipy.integrate import simps
+                        from scipy.integrate import simpson
 
                         axis = ax.axis[:ind]
-                        return simps(y=data, x=axis)
+                        return simpson(y=data, x=axis)
 
             I0 = self.map(
                 estimating_function,
@@ -1058,7 +1055,7 @@ def fourier_ratio_deconvolution(
         axis = ll.axes_manager.signal_axes[0]
         if fwhm is None:
             fwhm = float(
-                ll.get_current_signal().estimate_peak_width()._get_current_data()
+                ll.get_current_signal().estimate_peak_width()._get_current_data()[0]
             )
             _logger.info("FWHM = %1.2f" % fwhm)
 
@@ -1069,9 +1066,7 @@ def fourier_ratio_deconvolution(
             I0_shape.insert(axis.index_in_array, 1)
             I0 = I0.reshape(I0_shape)
 
-        from hyperspy.components1d import Gaussian
-
-        g = Gaussian()
+        g = hs.model.components1D.Gaussian()
         g.sigma.value = fwhm / 2.3548
         g.A.value = 1
         g.centre.value = 0
@@ -1287,7 +1282,7 @@ def power_law_extrapolation(
             s.data = np.zeros(new_shape)
             s.data[..., : axis.size] = self.data
         s.get_dimensions_from_data()
-        pl = PowerLaw()
+        pl = hs.model.components1D.PowerLaw()
         pl._axes_manager = self.axes_manager
         A, r = pl.estimate_parameters(
             s,
@@ -1454,7 +1449,7 @@ def kramers_kronig_analysis(
         axis = s.axes_manager.signal_axes[0]
         eaxis = axis.axis.copy()
 
-        if isinstance(zlp, hyperspy.signal.BaseSignal):
+        if isinstance(zlp, hs.signals.BaseSignal):
             if (
                 zlp.axes_manager.navigation_dimension
                 == self.axes_manager.navigation_dimension
@@ -1480,7 +1475,7 @@ def kramers_kronig_analysis(
                              in the BaseSignal class or a Number."
             )
 
-        if isinstance(t, hyperspy.signal.BaseSignal):
+        if isinstance(t, hs.signals.BaseSignal):
             if (
                 t.axes_manager.navigation_dimension
                 == self.axes_manager.navigation_dimension
@@ -1711,13 +1706,13 @@ def _get_offsets_and_segments(self, edges):
         return offsets, segments
 
     def _initialise_markers(self):
-        self._edge_markers["lines"] = Lines(
+        self._edge_markers["lines"] = hs.plot.markers.Lines(
             segments=np.empty((0, 2, 2)),
             transform="relative",
             color="black",
             shift=np.array([0.0, 0.19]),
         )
-        self._edge_markers["texts"] = Texts(
+        self._edge_markers["texts"] = hs.plot.markers.Texts(
             offsets=np.empty((0, 2)),
             texts=np.empty((0,)),
             offset_transform="relative",
@@ -1975,7 +1970,7 @@ def rebin(self, new_shape=None, scale=None, crop=True, dtype=None, out=None):
             out.events.data_changed.trigger(obj=out)
         return m
 
-    rebin.__doc__ = hyperspy.signal.BaseSignal.rebin.__doc__
+    rebin.__doc__ = hs.signals.BaseSignal.rebin.__doc__
 
     def vacuum_mask(
         self, threshold=10.0, start_energy=None, closing=True, opening=False
diff --git a/exspy/tests/models/test_linear_model.py b/exspy/tests/models/test_linear_model.py
index b81d8e951..d5f503702 100644
--- a/exspy/tests/models/test_linear_model.py
+++ b/exspy/tests/models/test_linear_model.py
@@ -19,8 +19,6 @@
 import numpy as np
 import pytest
 from hyperspy.decorators import lazifyTestClass
-from hyperspy._components.gaussian import Gaussian
-from hyperspy._components.lorentzian import Lorentzian
 import hyperspy.api as hs
 
 import exspy
@@ -100,8 +98,8 @@ def setup_method(self, method):
         s = hs.data.two_gaussians().inav[0]
         s.set_signal_type("EELS")
         m = s.create_model(auto_background=False, auto_add_edges=False)
-        g1 = Gaussian(centre=40)
-        g2 = Gaussian(centre=55)
+        g1 = hs.model.components1D.Gaussian(centre=40)
+        g2 = hs.model.components1D.Gaussian(centre=55)
         m.extend([g1, g2])
 
         # make dummy twinning
@@ -134,7 +132,7 @@ def test_expression_convolved(nav_dim, multiple_free_parameters):
     to_convolve.axes_manager[-1].offset = -to_convolve_component.centre.value
 
     # Create reference signal from model with convolution
-    l_ref = Lorentzian(A=100, centre=20, gamma=4)
+    l_ref = hs.model.components1D.Lorentzian(A=100, centre=20, gamma=4)
     m_ref = s_ref.create_model(auto_add_edges=False, auto_background=False)
     m_ref.append(l_ref)
     m_ref.low_loss = to_convolve
@@ -148,7 +146,7 @@ def test_expression_convolved(nav_dim, multiple_free_parameters):
         to_convolve = hs.stack([to_convolve] * 3)
 
     m = s.create_model(auto_add_edges=False, auto_background=False)
-    lor = Lorentzian(centre=20, gamma=4)
+    lor = hs.model.components1D.Lorentzian(centre=20, gamma=4)
     m.append(lor)
     assert not m.convolved
     m.low_loss = to_convolve
@@ -179,7 +177,7 @@ def test_expression_multiple_linear_parameter(nav_dim, convolve):
     p_ref = hs.model.components1D.Polynomial(order=2, a0=25, a1=-50, a2=2.5)
 
     # Create signal to convolve
-    to_convolve_component = Gaussian(A=100, sigma=5, centre=10)
+    to_convolve_component = hs.model.components1D.Gaussian(A=100, sigma=5, centre=10)
     to_convolve = hs.signals.Signal1D(to_convolve_component.function(np.arange(1000)))
     to_convolve.axes_manager[-1].offset = -to_convolve_component.centre.value
 
@@ -224,12 +222,12 @@ def test_multiple_linear_parameters_convolution(nav_dim):
     s_ref = EELSSpectrum(np.ones(1000))
 
     # Create signal to convolve
-    to_convolve_component = Gaussian(A=1000, sigma=50, centre=100)
+    to_convolve_component = hs.model.components1D.Gaussian(A=1000, sigma=50, centre=100)
     to_convolve = EELSSpectrum(to_convolve_component.function(np.arange(1000)))
     to_convolve.axes_manager[-1].offset = -to_convolve_component.centre.value
 
-    l_ref1 = Lorentzian(A=100, centre=200, gamma=10)
-    l_ref2 = Lorentzian(A=100, centre=600, gamma=20)
+    l_ref1 = hs.model.components1D.Lorentzian(A=100, centre=200, gamma=10)
+    l_ref2 = hs.model.components1D.Lorentzian(A=100, centre=600, gamma=20)
 
     m_ref = s_ref.create_model(auto_add_edges=False, auto_background=False)
     m_ref.extend([l_ref1, l_ref2])
@@ -244,8 +242,8 @@ def test_multiple_linear_parameters_convolution(nav_dim):
         to_convolve = hs.stack([to_convolve] * 3)
 
     m = s.create_model(auto_add_edges=False, auto_background=False)
-    l1 = Lorentzian(centre=200, gamma=10)
-    l2 = Lorentzian(centre=600, gamma=20)
+    l1 = hs.model.components1D.Lorentzian(centre=200, gamma=10)
+    l2 = hs.model.components1D.Lorentzian(centre=600, gamma=20)
     m.extend([l1, l2])
     assert not m.convolved
     m.low_loss = to_convolve
diff --git a/upcoming_changes/33.maintenance.rst b/upcoming_changes/33.maintenance.rst
new file mode 100644
index 000000000..8c3e1c7d7
--- /dev/null
+++ b/upcoming_changes/33.maintenance.rst
@@ -0,0 +1 @@
+Fix deprecation scipy and numpy warnings.
\ No newline at end of file