diff --git a/lcs/__init__.py b/lcs/__init__.py
index c20a3a1..73fbc26 100644
--- a/lcs/__init__.py
+++ b/lcs/__init__.py
@@ -1,4 +1,5 @@
-from . import contagion, inference, utilities
+from . import contagion, generative, inference, utilities
 from .contagion import *
+from .generative import *
 from .inference import *
 from .utilities import *
diff --git a/lcs/generative.py b/lcs/generative.py
index 3a5afbb..4e05e59 100644
--- a/lcs/generative.py
+++ b/lcs/generative.py
@@ -1,13 +1,28 @@
 import random
 
 import numpy as np
+import networkx as nx
+import xgi
 
 
-def erdos_renyi(n, p):
-    A = np.zeros((n, n))
+def erdos_renyi(n, p, seed=None):
+    return nx.adjacency_matrix(nx.fast_gnp_random_graph(n, p, seed)).todense()
 
-    for i in range(n):
-        for j in range(i):
-            A[i, j] = A[j, i] = random.random() <= p
 
-    return A
+def watts_strogatz(n, k, p, seed=None):
+    return nx.adjacency_matrix(nx.watts_strogatz_graph(n, k, p, seed)).todense()
+
+
+def sbm(n, k, epsilon, seed=None):
+    p = k / (n - 1)
+    # ratio of inter- to intra-community edges
+    p_in = (1 + epsilon) * p
+    p_out = (1 - epsilon) * p
+    return nx.adjacency_matrix(
+        nx.planted_partition_graph(2, int(n / 2), p_in, p_out, seed=seed)
+    ).todense()
+
+
+def projected_bipartite(k, s, seed=None):
+    H = xgi.chung_lu_hypergraph(k, s, seed)
+    return xgi.adjacency_matrix(H, sparse=False)
diff --git a/lcs/utilities.py b/lcs/utilities.py
index 611f886..2b0b753 100644
--- a/lcs/utilities.py
+++ b/lcs/utilities.py
@@ -37,8 +37,13 @@ def hamming_distance(A1, A2):
     return np.sum(np.abs(A1 - A2)) / 2
 
 
-def infections_per_node(x):
-    return np.mean(np.sum(x[1:] - x[:-1] > 0, axis=0))
+def infections_per_node(x, mode="mean"):
+    if mode == "mean":
+        return np.mean(np.sum(x[1:] - x[:-1] > 0, axis=0))
+    if mode == "median":
+        return np.median(np.sum(x[1:] - x[:-1] > 0, axis=0))
+    if mode == "max":
+        return np.max(np.sum(x[1:] - x[:-1] > 0, axis=0))
 
 
 def nu_distribution(x, A):
@@ -57,3 +62,17 @@ def degrees(A):
     if not isinstance(A, np.ndarray):
         A = A.todense()
     return A.sum(axis=0)
+
+
+def powerlaw(n, minval, maxval, r):
+    u = np.random.random(n)
+    return (minval**(1-r) + u*(maxval**(1-r) - minval**(1-r)))**(1/(1-r))
+
+
+def mean_power_law(minval, maxval, r):
+    if r == 1:
+        return -(minval - maxval) / (np.log(maxval) - np.log(minval))
+    elif r == 2:
+        return (np.log(maxval) - np.log(minval)) / (1/minval - 1/maxval)
+    else:
+        return (minval**(2-r)-maxval**(2-r))*(r-1)/((minval**(1-r)-maxval**(1-r))*(r-2))
\ No newline at end of file
diff --git a/requirements.txt b/requirements.txt
index 4becd88..f1ea0a2 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,6 +1,7 @@
 numpy
 scipy
-networkx
+networkx>=3.0
 matplotlib
 numba
-pytest
\ No newline at end of file
+pytest
+xgi
\ No newline at end of file