diff --git a/deerlab/bootstrap_analysis.py b/deerlab/bootstrap_analysis.py
index fdcbd3f35..d65e3fcf5 100644
--- a/deerlab/bootstrap_analysis.py
+++ b/deerlab/bootstrap_analysis.py
@@ -5,10 +5,9 @@
 
 import numpy as np
 import types
-from tqdm.auto import tqdm
-from joblib import Parallel, delayed
+from joblib import delayed
 from deerlab.classes import UQResult
-from deerlab.utils import isnumeric
+from deerlab.utils import isnumeric, _ProgressParallel
 from deerlab.noiselevel import noiselevel
 
 def bootstrap_analysis(fcn,Vexp,Vfit, samples=1000, noiselvl=None, resampling='gaussian', verbose = False, cores=1, memorylimit=8):
@@ -188,25 +187,5 @@ def bootsample():
 
 
 #-------------------------------------------------------------------------------
-class _ProgressParallel(Parallel):
-    """
-    Patch for joblib.Parallel
 
-    Overrides the print_progress() method to enable the synchronous use of the TQDM bar
-    even for parallel processing.  
-    """
-    def __init__(self, use_tqdm=True, total=None, *args, **kwargs):
-        self._use_tqdm = use_tqdm
-        self._total = total
-        super().__init__(*args, **kwargs)
-
-    def __call__(self, *args, **kwargs):
-        with tqdm(disable=not self._use_tqdm, total=self._total) as self._pbar:
-            return Parallel.__call__(self, *args, **kwargs)
-
-    def print_progress(self):
-        if self._total is None:
-            self._pbar.total = self.n_dispatched_tasks
-        self._pbar.n = self.n_completed_tasks
-        self._pbar.refresh()
 #-------------------------------------------------------------------------------
\ No newline at end of file
diff --git a/deerlab/profile_analysis.py b/deerlab/profile_analysis.py
index c67d6febb..8de93e2a4 100644
--- a/deerlab/profile_analysis.py
+++ b/deerlab/profile_analysis.py
@@ -7,10 +7,12 @@
 from scipy.stats import chi2 
 from deerlab import fit
 from deerlab import noiselevel,UQResult
+from deerlab.utils import _ProgressParallel
 import warnings 
 from tqdm import tqdm
+from joblib import delayed
 
-def profile_analysis(model,y, *args, parameters='all', grids=None, samples=50, noiselvl=None, verbose=False,**kargs):
+def profile_analysis(model,y, *args, parameters='all', grids=None, samples=50, noiselvl=None, verbose=False, cores=1, **kargs):
     r""" 
     Profile likelihood analysis for uncertainty quantification
 
@@ -49,6 +51,11 @@ def profile_analysis(model,y, *args, parameters='all', grids=None, samples=50, n
         Specifies whether to print the progress of the bootstrap analysis on the 
         command window, the default is false.
 
+    cores : scalar, optional
+        Number of CPU cores/processes for parallel computing. If ``cores=1`` no parallel 
+        computing is used. If ``cores=-1`` all available CPUs are used. The default is 
+        one core (no parallelization).
+
     kargs : keyword-argument pairs
         Any other keyword-argument pairs to be passed to the ``fit`` function. See the 
         documentation of the ``fit`` function for further details. 
@@ -106,9 +113,7 @@ def profile_analysis(model,y, *args, parameters='all', grids=None, samples=50, n
             tqdm.write(f"Profiling model parameter '{parameter}':",end='')
 
         # Calculate the profile objective function for the parameter
-        profile = np.zeros(len(grid))
-        for n,value in enumerate(tqdm(grid, disable=not verbose)): 
-
+        def profile_sample(value):
             # Freeze the model parameter at current value
             getattr(model, parameter).freeze(value)
 
@@ -116,12 +121,12 @@ def profile_analysis(model,y, *args, parameters='all', grids=None, samples=50, n
             with warnings.catch_warnings():
                 warnings.simplefilter("ignore")
                 fitresult_ = fit(model, y, *args, **kargs)
-
-            # Extract the objective function value
-            profile[n] = fitresult_.cost
-
+            cost = fitresult_.cost
             # Unfreeze the parameter
             getattr(model, parameter).unfreeze()
+            return cost
+        
+        profile = _ProgressParallel(n_jobs=cores,total=len(grid),use_tqdm=verbose)(delayed(profile_sample)(value) for value in grid)
 
         profile = {'x':np.squeeze(grid),'y':profile}
         uqresults[parameter] = UQResult('profile', data=getattr(fitresult,parameter), profiles=profile, threshold=threshold, noiselvl=noiselvl)
diff --git a/deerlab/utils.py b/deerlab/utils.py
index 46b139dce..89ff4f68f 100644
--- a/deerlab/utils.py
+++ b/deerlab/utils.py
@@ -6,6 +6,8 @@
 from types import FunctionType
 from functools import wraps
 import pickle
+from tqdm.auto import tqdm
+from joblib import Parallel
 
 def parse_multidatasets(y_, A, weights, noiselvl, precondition=False, masks=None, subsets=None):
 #===============================================================================
@@ -282,7 +284,6 @@ def goodness_of_fit(x,xfit,Ndof,noiselvl):
             stats['aic'] - Akaike information criterion
             stats['aicc'] - Corrected Akaike information criterion
             stats['bic'] - Bayesian information criterion
-            stats['autocorr'] - Autocorrelation based on Durbin–Watson statistic
     """
     sigma = noiselvl
     Ndof = np.maximum(Ndof,1)
@@ -956,6 +957,29 @@ def choleskycovmat(Q,cholfactors):
     return Σ 
 # ----------------------------------------------------------------------------------
 
+class _ProgressParallel(Parallel):
+    """
+    Patch for joblib.Parallel
+
+    Overrides the print_progress() method to enable the synchronous use of the TQDM bar
+    even for parallel processing.  
+    """
+    def __init__(self, use_tqdm=True, total=None, *args, **kwargs):
+        self._use_tqdm = use_tqdm
+        self._total = total
+        super().__init__(*args, **kwargs)
+
+    def __call__(self, *args, **kwargs):
+        with tqdm(disable=not self._use_tqdm, total=self._total) as self._pbar:
+            return Parallel.__call__(self, *args, **kwargs)
+
+    def print_progress(self):
+        if self._total is None:
+            self._pbar.total = self.n_dispatched_tasks
+        self._pbar.n = self.n_completed_tasks
+        self._pbar.refresh()
+        
+# ----------------------------------------------------------------------------------
 
 def _config():
     import importlib 
diff --git a/examples/advanced/ex_profileanalysis.py b/examples/advanced/ex_profileanalysis.py
index 6a195bd50..def941c26 100644
--- a/examples/advanced/ex_profileanalysis.py
+++ b/examples/advanced/ex_profileanalysis.py
@@ -4,6 +4,12 @@
 -------------------------------------------------------------------
 
 How to obtain objective function profiles for the non-linear parameters of a model. 
+
+Profile analysis works by running the fitting routine multiple times while freezing the parameter of interest at different values. 
+This way, the profile of the objective function for this parameter can be obtained. 
+The confidence interval can then be estimated by comparing the profile to a statistical threshold. 
+To keep analysis time resonable, the number of samples for the profile can be specified via the ``samples`` keyword. 
+Additionally, the number of CPU cores for parallelization can be specified via the ``cores`` keyword.
 """ 
 
 import numpy as np
@@ -43,7 +49,7 @@
 print(results)
 
 # Compute uncertainty with the likelihood profile method for the spin concentration and modulation depth parameters
-profile_uq = dl.profile_analysis(Vmodel,Vexp,samples=20, parameters=['conc','mod']) 
+profile_uq = dl.profile_analysis(Vmodel,Vexp,samples=20, parameters=['conc','mod'], cores=4) 
 
 #%%
 # Extract fitted dipolar signal
diff --git a/examples/basic/ex_bootstrapping.py b/examples/basic/ex_bootstrapping.py
index 196bd772e..1e227ad1c 100644
--- a/examples/basic/ex_bootstrapping.py
+++ b/examples/basic/ex_bootstrapping.py
@@ -13,7 +13,9 @@
 confidence intervals to get accurate estimates of the uncertainty. 
 Conviniently, ``fit`` integrates bootstrapping to make it accessible via the keyword argument ``bootstrap`` which specifies the
 number of samples to analyze to estimate the uncertainty. The larger this number, the more accurate 
-the confidence intervals but the longer the analysis will be. The standard for publication is typically 1000 samples. 
+the confidence intervals but the longer the analysis will be. The standard for publication is typically 250 samples. 
+
+To keep the analysis time reasonable, bootstrapping can be parallelized by specifying the number of CPU cores to use via the ``bootcores`` keyword.
 """ 
 
 import numpy as np
@@ -48,7 +50,7 @@
 Vmodel = dl.dipolarmodel(t,r, experiment = dl.ex_4pdeer(tau1,tau2, pathways=[1]))
 
 # Fit the model to the data
-results = dl.fit(Vmodel,Vexp,bootstrap=20)
+results = dl.fit(Vmodel,Vexp,bootstrap=20,bootcores=4)
 
 # In this example, just for the sake of time, we will just use 20 bootstrap samples.  
 
diff --git a/test/test_profile_analysis.py b/test/test_profile_analysis.py
index 6be29fdc4..abe54b772 100644
--- a/test/test_profile_analysis.py
+++ b/test/test_profile_analysis.py
@@ -76,6 +76,25 @@ def test_grids(model, mock_data):
     assert np.allclose(x,grid['mean'])
 # ======================================================================
 
+# ======================================================================
+def test_ncores(model, mock_data):
+    "Check that the number of cores can be specified"
+    from unittest.mock import patch
+    from deerlab.utils import _ProgressParallel
+
+    n_jobs_used = []
+
+    class TrackingParallel(_ProgressParallel):
+        def __init__(self, *args, **kwargs):
+            n_jobs_used.append(kwargs.get('n_jobs'))
+            super().__init__(*args, **kwargs)
+
+    with patch('deerlab.profile_analysis._ProgressParallel', TrackingParallel):
+        profile_analysis(model, mock_data, samples=5, noiselvl=noiselvl, cores=4)
+
+    assert any(n == 4 for n in n_jobs_used)
+# ======================================================================
+
 # ======================================================================
 def test_docstring():
     "Check that the docstring includes all variables and keywords."