diff --git a/tests/testthat/test-feature-unphased.R b/tests/testthat/test-feature-unphased.R
index 024805a..43aba79 100644
--- a/tests/testthat/test-feature-unphased.R
+++ b/tests/testthat/test-feature-unphased.R
@@ -39,17 +39,38 @@ test_that("generating the scrm command works", {
 
 test_that("unphasing works", {
   seg_sites <- list()
-  seg_sites[[1]] <-   matrix(c(1, 0, 0, 0,
-                               1, 1, 0, 1,
-                               1, 0, 0, 1,
-                               1, 0, 0, 1), 4, 4, byrow = TRUE)
-  attr(seg_sites[[1]], "positions") <- c(0.1, 0.2, 0.5, 0.7)
-  attr(seg_sites[[1]], "locus") <- rep(0, each = 4)
+  seg_sites[[1]] <- create_segsites(matrix(c(0, 1, 0, 1,
+                                             1, 0, 1, 0,
+                                             1, 0, 1, 1,
+                                             0, 1, 0, 0), 4, 4, byrow = TRUE),
+                                    c(0.1, 0.2, 0.5, 0.7))
+
 
   phased <- unphase_segsites(seg_sites, 2, 1)
   expect_that(phased, is_a("list"))
   expect_that(dim(phased[[1]])[1], is_equivalent_to(2))
-  expect_that(dim(phased[[1]])[2], is_less_than(3))
+  expect_that(dim(phased[[1]])[2], is_less_than(5))
+
+  n_snps <- sapply(1:10000, function(i) {
+    phased <- unphase_segsites(seg_sites, 2, 1)
+    c(ncol(phased[[1]]), sum(phased[[1]][1, ]))
+  })
+
+  expect_less_than(sum(abs(table(n_snps[1 , ]) / ncol(n_snps) -
+                             dbinom(0:4, 4, .5))), 0.1)
+  expect_less_than(sum(abs(table(n_snps[2 , ]) / ncol(n_snps) -
+                             dbinom(0:4, 4, .25))), 0.1)
+
+
+  seg_sites[[1]] <- create_segsites(matrix(c(0, 1, 0, 1,
+                                             1, 0, 1, 0,
+                                             1, 1, 1, 1,
+                                             1, 1, 1, 1), 4, 4, byrow = TRUE),
+                                    c(0.1, 0.2, 0.5, 0.7))
+  phased <- unphase_segsites(seg_sites, 2, 1)
+  expect_true(all(phased[[1]][1, ] == 0))
+  expect_true(all(phased[[1]][2, ] == 1))
+
 
   phased <- unphase_segsites(seg_sites, 2, 2)
   expect_that(phased, is_a("list"))
@@ -57,10 +78,6 @@ test_that("unphasing works", {
   expect_equal(dim(phased[[1]]), c(4, 4))
   expect_equal(colSums(seg_sites[[1]]), colSums(phased[[1]]))
 
-  expect_equal(attr(phased[[1]], "positions"),
-               attr(seg_sites[[1]], "positions"))
-  expect_equal(attr(phased[[1]], "locus"),
-               attr(seg_sites[[1]], "locus"))
 
   seg_sites[[2]] <- seg_sites[[1]]
   phased <- unphase_segsites(seg_sites, 2, 1)
@@ -69,9 +86,17 @@ test_that("unphasing works", {
 
   seg_sites[[1]] <- seg_sites[[1]][ , numeric()]
   attr(seg_sites[[1]], "positions") <- numeric()
-  attr(seg_sites[[1]], "locus") <- numeric()
+  #attr(seg_sites[[1]], "locus") <- numeric()
   phased <- unphase_segsites(seg_sites, 2, 1)
-  expect_equivalent(seg_sites[[1]], matrix(0, 4, 0))
+  #expect_equivalent(seg_sites[[1]], matrix(0, 4, 0))
   expect_equal(attr(seg_sites[[1]], "positions"), numeric(0))
-  expect_equal(attr(seg_sites[[1]], "locus"), numeric(0))
+  #expect_equal(attr(seg_sites[[1]], "locus"), numeric(0))
+})
+
+
+test_that("simulating unphased data works", {
+  scrm <- get_simulator("scrm")
+  model <- coal_model(5, 1) + feat_unphased(2, 1) + feat_mutation(5) + sumstat_seg_sites()
+  data <- simulate(model)
+  expect_equal(nrow(data$seg_sites[[1]]), 5)
 })