diff --git a/hera_pspec/pstokes.py b/hera_pspec/pstokes.py
index fd3d33b9..8a0d33fd 100644
--- a/hera_pspec/pstokes.py
+++ b/hera_pspec/pstokes.py
@@ -188,7 +188,7 @@ def _combine_pol_arrays(
 
     nsamples_list :  any length 2 iterable of numpy arrays
         Iterable of nsamples arrays to be combined to form their pseudo-Stokes equivalent.
-        nsamples are combined to preserve proper variance, see hera_pspec issue #391.
+        nsamples are combined to preserve proper variance, see hera_pspec issue #391 and #406
         If only one is given, it is duplicated.
         Default is None.
 
@@ -264,9 +264,9 @@ def _combine_pol_arrays(
 
     # constructing nsamples
     if nsamples_list is not None:
-        combined_nsamples = (nsamples_list[0]**-1 + nsamples_list[1]**-1)**-1
-        if pol_convention == 'avg':
-            combined_nsamples *= 4.
+        combined_nsamples = 4 * (nsamples_list[0]**-1 + nsamples_list[1]**-1)**-1
+        # for an explanation of the factor of 4, and why it doesn't depend on pol_convention,
+        # see https://github.com/HERA-Team/hera_pspec/issues/406
     else:
         combined_nsamples = None
 
diff --git a/hera_pspec/tests/test_pstokes.py b/hera_pspec/tests/test_pstokes.py
index d5553389..fd870598 100644
--- a/hera_pspec/tests/test_pstokes.py
+++ b/hera_pspec/tests/test_pstokes.py
@@ -20,8 +20,10 @@ def setUp(self):
         # Loading pyuvdata objects
         self.uvd1 = pyuvdata.UVData()
         self.uvd1.read_miriad(dset1, use_future_array_shapes=True)
+        setattr(self.uvd1, 'vis_units', 'Jy')
         self.uvd2 = pyuvdata.UVData()
         self.uvd2.read_miriad(dset2, use_future_array_shapes=True)
+        setattr(self.uvd2, 'vis_units', 'Jy')
 
     def tearDown(self):
         pass
@@ -43,9 +45,10 @@ def test_combine_pol(self):
         # out1 assumed avg by default
         setattr(uvd1, 'pol_convention', 'sum')
         setattr(uvd2, 'pol_convention', 'sum')
+        
         out3 = pstokes._combine_pol(uvd1, uvd2, 'XX', 'YY')
         assert np.allclose(out3.data_array, out1.data_array * 2.)
-        assert np.allclose(out3.nsample_array, out1.nsample_array / 4.)
+        assert np.allclose(out3.nsample_array, out1.nsample_array)
 
         # check exceptions
         pytest.raises(AssertionError, pstokes._combine_pol, dset1, dset2, 'XX', 'YY' )
@@ -81,7 +84,7 @@ def test_combine_pol_arrays(self):
         assert d2 is None
         assert f2 is None
         assert ns2 is None
-        # polarizations can be strings     
+        # polarizations can be strings
         d3, f3, ns3 = pstokes._combine_pol_arrays(
             pol1='XX',
             pol2='YY',
@@ -153,6 +156,7 @@ def test_construct_pstokes_multipol(self):
         uvd.read(multipol_dset, use_future_array_shapes=True)
         uvc = pyuvdata.UVCal()
         uvc.read_calfits(multipol_dset_cal, use_future_array_shapes=True)
+        uvc.gain_scale = 'Jy'
         uvutils.uvcalibrate(uvd, uvc)
         wgts = [(0.5, 0.5), (0.5, -0.5)]