diff --git a/examples/injection_examples/basic_tutorial.py b/examples/injection_examples/basic_tutorial.py
index f9946da0039934f222a9072026a2f773c67d0ea1..0d5ede77386bc38853aba7ca63042967a5d51efb 100644
--- a/examples/injection_examples/basic_tutorial.py
+++ b/examples/injection_examples/basic_tutorial.py
@@ -29,10 +29,10 @@ injection_parameters = dict(mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5
                             waveform_approximant='IMRPhenomPv2', reference_frequency=50., ra=1.375, dec=-1.2108)
 
 # Create the waveform_generator using a LAL BinaryBlackHole source function
-waveform_generator = tupak.waveform_generator.WaveformGenerator(time_duration=time_duration,
-                                                                sampling_frequency=sampling_frequency,
-                                                                frequency_domain_source_model=tupak.source.lal_binary_black_hole,
-                                                                parameters=injection_parameters)
+waveform_generator = tupak.WaveformGenerator(time_duration=time_duration,
+                                             sampling_frequency=sampling_frequency,
+                                             frequency_domain_source_model=tupak.source.lal_binary_black_hole,
+                                             parameters=injection_parameters)
 hf_signal = waveform_generator.frequency_domain_strain()
 
 # Set up interferometers.  In this case we'll use three interferometers (LIGO-Hanford (H1), LIGO-Livingston (L1),
@@ -55,11 +55,11 @@ for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'phase', 'psi'
 priors['luminosity_distance'] = tupak.prior.create_default_prior(name='luminosity_distance')
 
 # Initialise the likelihood by passing in the interferometer data (IFOs) and the waveoform generator
-likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator)
+likelihood = tupak.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator)
 
 # Run sampler.  In this case we're going to use the `dynesty` sampler
-result = tupak.sampler.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,
-                                   injection_parameters=injection_parameters, outdir=outdir, label=label)
+result = tupak.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,
+                           injection_parameters=injection_parameters, outdir=outdir, label=label)
 
 # make some plots of the outputs
 result.plot_corner()
diff --git a/examples/open_data_examples/GW150914.py b/examples/open_data_examples/GW150914.py
index 56ea1cde47d0274b68c2088ccd79ca571937d5a2..49eb90fabb2277e21be44da9344011bf2110bba2 100644
--- a/examples/open_data_examples/GW150914.py
+++ b/examples/open_data_examples/GW150914.py
@@ -44,20 +44,19 @@ prior['luminosity_distance'] = tupak.prior.PowerLaw(
 # creates the frequency-domain strain. In this instance, we are using the
 # `lal_binary_black_hole model` source model. We also pass other parameters:
 # the waveform approximant and reference frequency.
-waveform_generator = tupak.waveform_generator.WaveformGenerator(time_duration=interferometers[0].duration,
-                                                                sampling_frequency=interferometers[
-                                                                    0].sampling_frequency,
-                                                                frequency_domain_source_model=tupak.source.lal_binary_black_hole,
-                                                                parameters={'waveform_approximant': 'IMRPhenomPv2',
-                                                                            'reference_frequency': 50})
+waveform_generator = tupak.WaveformGenerator(time_duration=interferometers[0].duration,
+                                             sampling_frequency=interferometers[0].sampling_frequency,
+                                             frequency_domain_source_model=tupak.source.lal_binary_black_hole,
+                                             parameters={'waveform_approximant': 'IMRPhenomPv2',
+                                                         'reference_frequency': 50})
 
 # In this step, we define the likelihood. Here we use the standard likelihood
 # function, passing it the data and the waveform generator.
-likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers, waveform_generator)
+likelihood = tupak.GravitationalWaveTransient(interferometers, waveform_generator)
 
 # Finally, we run the sampler. This function takes the likelihood and prio
 # along with some options for how to do the sampling and how to save the data
-result = tupak.sampler.run_sampler(likelihood, prior, sampler='dynesty',
-                                   outdir=outdir, label=label)
+result = tupak.run_sampler(likelihood, prior, sampler='dynesty',
+                           outdir=outdir, label=label)
 result.plot_corner()
 print(result)
diff --git a/examples/open_data_examples/GW150914_minimal.py b/examples/open_data_examples/GW150914_minimal.py
index 606409387b622ffc583dff647af6bebd94b01603..7df645cbbd093c3856920a67cbbbbe05e632dd05 100644
--- a/examples/open_data_examples/GW150914_minimal.py
+++ b/examples/open_data_examples/GW150914_minimal.py
@@ -10,5 +10,5 @@ t0 = tupak.utils.get_event_time("GW150914")
 prior = dict(geocent_time=tupak.prior.Uniform(t0-0.1, t0+0.1, name='geocent_time'))
 interferometers = tupak.detector.get_event_data("GW150914")
 likelihood = tupak.likelihood.get_binary_black_hole_likelihood(interferometers)
-result = tupak.sampler.run_sampler(likelihood, prior, label='GW150914')
+result = tupak.run_sampler(likelihood, prior, label='GW150914')
 result.plot_corner()
diff --git a/examples/other_examples/linear_regression.py b/examples/other_examples/linear_regression.py
index 1aee5f0243bdee5180bc384c09aa70cfe83a950e..208e05b97e63f7168ccf5b30df7f79af7cde5105 100644
--- a/examples/other_examples/linear_regression.py
+++ b/examples/other_examples/linear_regression.py
@@ -55,7 +55,7 @@ fig.savefig('{}/{}_data.png'.format(outdir, label))
 # our model.
 
 
-class GaussianLikelihood(tupak.likelihood.Likelihood):
+class GaussianLikelihood(tupak.Likelihood):
     def __init__(self, x, y, sigma, waveform_generator):
         """
 
@@ -91,7 +91,7 @@ class GaussianLikelihood(tupak.likelihood.Likelihood):
 # can generate a signal. We give it information on how to make the time series
 # and the model() we wrote earlier.
 
-waveform_generator = tupak.waveform_generator.WaveformGenerator(
+waveform_generator = tupak.WaveformGenerator(
     time_duration=time_duration, sampling_frequency=sampling_frequency,
     time_domain_source_model=model)
 
@@ -106,7 +106,7 @@ priors['m'] = tupak.prior.Uniform(0, 5, 'm')
 priors['c'] = tupak.prior.Uniform(-2, 2, 'c')
 
 # And run sampler
-result = tupak.sampler.run_sampler(
+result = tupak.run_sampler(
     likelihood=likelihood, priors=priors, sampler='dynesty', npoints=500,
     walks=10, injection_parameters=injection_parameters, outdir=outdir,
     label=label, plot=True)
diff --git a/tupak/__init__.py b/tupak/__init__.py
index 5716af07db5ea5e398f25e36ed030e502905e1f9..26003e316a429a19b505dbbadcfd6d5d7393920b 100644
--- a/tupak/__init__.py
+++ b/tupak/__init__.py
@@ -21,3 +21,8 @@ from . import waveform_generator
 from . import result
 from . import sampler
 from . import conversion
+
+# import a few oft-used functions and classes to simplify scripts
+from likelihood import Likelihood, GravitationalWaveTransient
+from waveform_generator import WaveformGenerator
+from sampler import run_sampler