diff --git a/distarray/dist/__init__.py b/distarray/dist/__init__.py
index 962fb631..fe04183e 100644
--- a/distarray/dist/__init__.py
+++ b/distarray/dist/__init__.py
@@ -9,4 +9,5 @@
 from distarray.dist.distarray import DistArray
 from distarray.dist.context import Context
 from distarray.dist.maps import Distribution
+from distarray.dist import random
 from distarray.dist.functions import *
diff --git a/distarray/fuzz/__init__.py b/distarray/fuzz/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/distarray/fuzz/_test_fuzz.py b/distarray/fuzz/_test_fuzz.py
new file mode 100755
index 00000000..87b9407d
--- /dev/null
+++ b/distarray/fuzz/_test_fuzz.py
@@ -0,0 +1,9 @@
+#!/usr/bin/env python
+
+import distarray.dist as da
+from distarray.fuzz.make_rand_distarray import rand_distarray
+from distarray.fuzz.fuzzy_test import fuzzy_test
+
+c = da.Context()
+for _ in range(2):
+    fuzzy_test(rand_distarray(c))
diff --git a/distarray/fuzz/fuzzy_test.py b/distarray/fuzz/fuzzy_test.py
new file mode 100644
index 00000000..485b55a9
--- /dev/null
+++ b/distarray/fuzz/fuzzy_test.py
@@ -0,0 +1,59 @@
+import random
+
+import numpy as np
+from numpy.testing import (assert_almost_equal, assert_array_almost_equal)
+
+import distarray.dist as dist
+from distarray.dist.decorators import local, vectorize
+
+
+def fuzzy_test(darr):
+    ndarr = darr.toarray()
+
+    def assert_dist_numpy_equal(darr, ndarr):
+        if darr is None and ndarr is None:
+            return
+        elif isinstance(darr, dist.DistArray):
+            assert_array_almost_equal(darr.toarray(), ndarr)
+        assert_almost_equal(darr, ndarr)
+
+    def test_method(darr, ndarr, meth_name, args=None, kwargs=None):
+        args = () if args is None else args
+        kwargs = {} if kwargs is None else kwargs
+        dist_meth = getattr(darr, meth_name)
+        nd_meth = getattr(ndarr, meth_name)
+        assert_dist_numpy_equal(dist_meth(*args, **kwargs),
+                                nd_meth(*args, **kwargs))
+
+    def random_index(darr):
+        shape = darr.shape
+        index = []
+        for i in range(len(shape)):
+            index.append(random.choice(range(shape[i])))
+        return tuple(index)
+
+    # test some reductions
+    test_method(darr, ndarr, 'sum')
+    test_method(darr, ndarr, 'var')
+    test_method(darr, ndarr, 'mean')
+    # test getting and setting
+    test_method(darr, ndarr, '__getitem__', args=(random_index(darr),))
+    test_method(darr, ndarr, '__setitem__', args=(random_index(darr),
+                                                  random.random()))
+    # test a binary function
+    assert_array_almost_equal(dist.add(darr, darr).toarray(), np.add(ndarr,
+                                                                     ndarr))
+    # test a unary function
+    assert_array_almost_equal(dist.cos(darr).toarray(), np.cos(ndarr))
+
+    # test local
+    @local
+    def local_add50(darr):
+        return darr + 50
+    assert_array_almost_equal(local_add50(darr).toarray(), ndarr + 50)
+
+    # test vectorize
+    @vectorize
+    def sqr_add10(a):
+        return a**2 + 10
+    assert_array_almost_equal(sqr_add10(darr).toarray(), ndarr**2 + 10)
diff --git a/distarray/fuzz/rand_distarray.py b/distarray/fuzz/rand_distarray.py
new file mode 100644
index 00000000..6b3cf17d
--- /dev/null
+++ b/distarray/fuzz/rand_distarray.py
@@ -0,0 +1,48 @@
+"""
+Generate a random distarray with reasonable defaults.
+"""
+
+import random
+
+import distarray.dist as da
+from distarray.dist import Distribution
+
+
+def rand_size(max_size=None):
+    if max_size is None:
+        max_size = 1e3
+    return random.randint(1, max_size)
+
+
+def rand_ndim(max_ndim=None):
+    if max_ndim is None:
+        max_ndim = 5
+    return random.randint(1, max_ndim)
+
+
+def rand_shape(max_1d=None, ndim=None):
+    if ndim is None:
+        ndim = rand_ndim()
+
+    if max_1d is None:
+        max_1d = 15
+    shape = []
+    for i in range(ndim):
+        shape.append(random.randint(1, max_1d))
+    return tuple(shape)
+
+
+def rand_dist(ndim):
+    choices = ('b', 'n', 'c')
+    dist = []
+    for i in range(ndim):
+        dist.append(random.choice(choices))
+    return tuple(dist)
+
+
+def rand_distarray(context):
+    shape = rand_shape()
+    dist = rand_dist(len(shape))
+    db = Distribution.from_shape(context, shape, dist=dist)
+    r = da.random.Random(context)
+    return r.randint(db, -1e3, 1e3)