diff --git a/hera_pspec/pstokes.py b/hera_pspec/pstokes.py
index 5cce59d5..f8fad5eb 100644
--- a/hera_pspec/pstokes.py
+++ b/hera_pspec/pstokes.py
@@ -145,7 +145,8 @@ def _combine_pol(uvd1, uvd2, pol1, pol2, pstokes='pI', x_orientation=None):
     uvdS.data_array = stdata  # pseudo-stokes data
     uvdS.flag_array = flag  # flag array
     uvdS.polarization_array = np.array([pstokes], dtype=int) # polarization number
-    uvdS.nsample_array = uvd1.nsample_array + uvd2.nsample_array # nsamples
+    # nsamples combined to preserve proper variance, see hera_pspec issue #391
+    uvdS.nsample_array = 4 * (uvd1.nsample_array**-1 + uvd2.nsample_array**-1)**-1 
 
     uvdS.history = "Merged into pseudo-stokes vis with hera_pspec version {}\n{}" \
                     "{}{}{}{}\n".format(__version__, "-"*20+'\n',