diff --git a/snipar/simulate.py b/snipar/simulate.py
index d4572221..eb4b0d72 100644
--- a/snipar/simulate.py
+++ b/snipar/simulate.py
@@ -470,7 +470,7 @@ def forward_sim(haps, maps, ngen_random, ngen_am, unlinked, n_causal, h2, v_indi
                 delta_p, delta_m, Y_p, Y_m = compute_phenotype(haps, causal, a, 1 - h2)
                 delta_p, delta_m, Y_p, Y_m = delta_p[father_indices], delta_m[mother_indices], Y_p[father_indices], Y_m[mother_indices]
             else:
-                delta_p, delta_m, Y_p, Y_m = G_males[father_indices], G_females[mother_indices], Y_males[father_indices], Y_females[father_indices]
+                delta_p, delta_m, Y_p, Y_m = G_males[father_indices], G_females[mother_indices], Y_males[father_indices], Y_females[mother_indices]
             # Compute phenotypes
             print('Computing phenotypes')
             if np.abs(beta_vert) > 0 and a_count>0:
@@ -541,7 +541,7 @@ def forward_sim(haps, maps, ngen_random, ngen_am, unlinked, n_causal, h2, v_indi
             delta_p, delta_m, Y_p, Y_m = compute_phenotype(haps, causal, a, 1 - h2)
             delta_p, delta_m, Y_p, Y_m = delta_p[father_indices], delta_m[mother_indices], Y_p[father_indices], Y_m[mother_indices]
         else:
-            delta_p, delta_m, Y_p, Y_m = G_males[father_indices], G_females[mother_indices], Y_males[father_indices], Y_females[father_indices]
+            delta_p, delta_m, Y_p, Y_m = G_males[father_indices], G_females[mother_indices], Y_males[father_indices], Y_females[mother_indices]
         # Compute phenotypes
         print('Computing phenotypes')
         if np.abs(beta_vert) > 0 and a_count>0:
@@ -629,4 +629,4 @@ def impute_all_fams_phased(haps,freqs,ibd):
     for i in prange(haps.shape[0]):
         for j in range(haps.shape[2]):
             imp[i,j] = impute_from_sibs_phased(haps[i,0,j,:],haps[i,1,j,:],ibd[i,j,:],freqs[j])
-    return imp
\ No newline at end of file
+    return imp