diff --git a/R/samples.R b/R/samples.R
new file mode 100644
index 00000000..ef3cb3d8
--- /dev/null
+++ b/R/samples.R
@@ -0,0 +1,49 @@
+# Three simple functions
+DEFAULT_N = 10000
+normal <- function (mean, std, n=DEFAULT_N){
+  return(rnorm(n, mean, std))
+}
+
+lognormal <- function(meanlog, sdlog, n=DEFAULT_N){
+  return(rlnorm(n, meanlog = meanlog, sdlog = sdlog))
+}
+
+to <- function(low, high, n=DEFAULT_N){
+  normal95confidencePoint = 1.6448536269514722
+  logLow = log(low)
+  logHigh = log(high)
+  meanlog = (logLow + logHigh)/2
+  sdlog = (logHigh - logLow) / (2.0 * normal95confidencePoint)
+  return(lognormal(meanlog, sdlog, n))
+}
+
+mixture <- function(samples_list, weights_array, n=DEFAULT_N){ # note that this takes a list, not an array
+  normalized_weights = weights_array/sum(weights_array)
+  cummulative_sums = cumsum(normalized_weights)
+  helper_probs = runif(n)
+  results = vector(mode='numeric', length=n)
+  for(i in c(1:n)){
+    helper_which_list = which(cummulative_sums > helper_probs[i])
+    helper_loc = ifelse(is.na(helper_which_list[1]), 1, helper_which_list[1])
+    target_samples = samples_list[[helper_loc]]
+    result = sample(target_samples, 1)
+    results[i] = result
+  }
+  return(results)
+}
+
+# Example
+p_a = 0.8
+p_b = 0.5
+p_c = p_a * p_b
+
+dists = list(c(0), c(1), to(1, 3), to(2, 10))
+# print(dists)
+weights = c((1 - p_c), p_c/2, p_c/4, p_c/4)
+# print(weights)
+result = mixture(dists, weights)
+mean_result = mean(result)
+print(mean_result)
+
+
+
diff --git a/README.md b/README.md
new file mode 100644
index 00000000..fe38f720
--- /dev/null
+++ b/README.md
@@ -0,0 +1 @@
+<https://github.com/rethinkpriorities/squigglepy>
\ No newline at end of file
diff --git a/example.squiggle b/example.squiggle
new file mode 100644
index 00000000..d83c6c3d
--- /dev/null
+++ b/example.squiggle
@@ -0,0 +1,6 @@
+p_a = 0.8
+p_b = 0.5
+p_c = p_a * p_b
+
+result = mx([0, 1, 1 to 3, 2 to 10], [(1 - p_c), p_c/2, p_c/4, p_c/4 ])
+mean(result)
\ No newline at end of file
diff --git a/python/samples.py b/python/samples.py
new file mode 100644
index 00000000..9e53d2fa
--- /dev/null
+++ b/python/samples.py
@@ -0,0 +1,55 @@
+# imports
+import numpy as np
+rng = np.random.default_rng(123)
+DEFAULT_N = 100000
+
+# three simple functions
+
+
+def normal(mean, std, n=DEFAULT_N):
+    return np.array(rng.normal(mean, std, n))
+
+
+def lognormal(mean, std, n=DEFAULT_N):
+    return np.array(rng.lognormal(mean, std, n))
+
+
+def to(low, high, n=DEFAULT_N):
+    normal95confidencePoint = 1.6448536269514722
+    logLow = np.log(low)
+    logHigh = np.log(high)
+    meanlog = (logLow + logHigh)/2
+    sdlog = (logHigh - logLow) / (2.0 * normal95confidencePoint)
+    return lognormal(meanlog, sdlog, n)
+
+
+def mixture(samples_list, weights_array, n=DEFAULT_N):
+    normalized_weights = weights_array/sum(weights_array)
+    cummulative_sums = np.cumsum(normalized_weights)
+    helper_probs = rng.random(n)
+    results = np.empty(n)
+    for i in range(n):
+        helper_list = [j for j in range(
+            len(cummulative_sums)) if cummulative_sums[j] > helper_probs[i]]
+        if len(helper_list) == 0:
+            helper_loc = 0
+        else:
+            helper_loc = helper_list[0]
+        target_samples = samples_list[helper_loc]
+        result = rng.choice(target_samples, 1)[0]
+        results[i] = result
+    return (results)
+
+
+# Example
+p_a = 0.8
+p_b = 0.5
+p_c = p_a * p_b
+
+dists = [[0], [1], to(1, 3), to(2, 10)]
+# print(dists)
+weights = np.array([1 - p_c, p_c/2, p_c/4, p_c/4])
+# print(weights)
+result = mixture(dists, weights)
+mean_result = np.mean(result)
+print(mean_result)