diff --git a/NAMESPACE b/NAMESPACE
index e44efa2..12ff81a 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -32,6 +32,7 @@ export(discount_function)
 export(invariance_checks)
 export(kirby_consistency)
 export(kirby_score)
+export(most_consistent_indiffs)
 export(nonsys)
 export(td_bclm)
 export(td_bcnm)
@@ -51,6 +52,7 @@ importFrom(stats,BIC)
 importFrom(stats,approx)
 importFrom(stats,binomial)
 importFrom(stats,coef)
+importFrom(stats,filter)
 importFrom(stats,fitted)
 importFrom(stats,glm)
 importFrom(stats,integrate)
diff --git a/R/generics.R b/R/generics.R
index 12d93f6..fe4220d 100644
--- a/R/generics.R
+++ b/R/generics.R
@@ -490,6 +490,8 @@ deviance.td_ddm <- function(mod) return(-2*logLik.td_ddm(mod))
 #' @param del Plots data for a particular delay
 #' @param val_del Plots data for a particular delayed value
 #' @param legend Logical: display a legend? Only relevant for \code{type = 'summary'} and \code{type = 'rt'}
+#' @param p_lines Numerical vector. When \code{type = 'summary'} the discount curve, where the probability of selecting the immediate reward is 0.5, is plotted. \code{p_lines} allows you to specify other probabilities for which similar curves should be plotted (only applicable for probabilistic models, e.g. \code{td_bcnm}, \code{td_bclm} and \code{td_ddm})
+#' @param p_tol If \code{p_lines} is not \code{NULL}, what is the maximum distance that estimated probabilities can be from their true values? Smaller values results in slower plot generation
 #' @param verbose Whether to print info about, e.g., setting del = ED50 when \code{type = 'endpoints'}
 #' @param confint When \code{type = 'rt'}, what confidence interval should be plotted for RTs? Default is 0.95 (95\% confidence interval)
 #' @param ... Additional arguments to \code{plot()}
@@ -497,13 +499,15 @@ deviance.td_ddm <- function(mod) return(-2*logLik.td_ddm(mod))
 #' \dontrun{
 #' data("td_bc_single_ptpt")
 #' mod <- td_bclm(td_bc_single_ptpt, model = 'hyperbolic.1')
-#' plot(mod, type = 'summary')
+#' plot(mod, type = 'summary', p_lines = c(0.25, 0.75), log = 'x')
 #' plot(mod, type = 'endpoints')
 #' }
 #' @export
 plot.td_um <- function(x,
                        type = c('summary', 'endpoints', 'link', 'rt'),
                        legend = T,
+                       p_lines = NULL,
+                       p_tol = 0.001,
                        verbose = T,
                        del = NULL,
                        val_del = NULL,
@@ -541,17 +545,42 @@ plot.td_um <- function(x,
     # Plot indifference curve
     lines(pred_indiffs ~ plotting_delays)
     
-    # Visualize stochasticity---goal for later. For now, don't know how to do this for td_bclm
-    # if (x$config$gamma_scale != 'none') {
-    #   if (verbose) {
-    #     cat(sprintf('gamma parameter (steepness of curve) is scaled by val_del.\nThus, the curve will have different steepness for a different value of val_del.\nDefaulting to val_del = %s (mean of val_del from data used to fit model).\nUse the `val_del` argument to specify a custom value.\n\n', val_del))
-    #   }
-    # }
-    # p_range <- args$p_range %def% c(0.4, 0.6)
-    # lower <- invert_decision_function(x, prob = p_range[1], del = plotting_delays)
-    # upper <- invert_decision_function(x, prob = p_range[2], del = plotting_delays)
-    # lines(lower ~ plotting_delays, lty = 'dashed')
-    # lines(upper ~ plotting_delays, lty = 'dashed')
+    if (!is.null(p_lines)) {
+      classes <- c('td_bcnm', 'td_bclm', 'td_ddm')
+      if (!inherits(x, classes)) {
+        stop(sprintf('p_lines can only be used for models of the following classes:\n%s',
+                     paste('-', classes, collapse = '\n')))
+      }
+      # Plot curves for other probabilities
+      
+      # Because of the overhead of individual calls to predict(), exhaustive grid search
+      # is faster than uniroot() or similar to invert predict()
+      
+      # Call predict() once on a big grid
+      val_imm_cands <- seq(0, val_del, length.out = ceiling(1/p_tol + 1))
+      grid <- data.frame(val_del = val_del,
+                         del = rep(plotting_delays, each = length(val_imm_cands)),
+                         val_imm = rep(val_imm_cands, times = length(plotting_delays)))
+      grid$p <- predict(x, grid, type = 'response')
+      
+      # Split the grid by delay
+      # Using split() with a numerical index is faster than calling tapply() or similar
+      split_idx <- rep(1:length(plotting_delays), each = length(val_imm_cands))
+      subgrid_list <- split(grid, split_idx)
+      for (p in p_lines) {
+        # Get the val_imm producing (close to) the desired p at each delay
+        val_imm <- sapply(subgrid_list, function(subgrid) {
+          if (max(subgrid$p) < p | min(subgrid$p) > p) {
+            return(NA)
+          } else {
+            return(subgrid$val_imm[which.min((subgrid$p - p)**2)])
+          }
+        })
+        # Plot
+        y <- val_imm / val_del
+        lines(y ~ plotting_delays, lty = 'dashed')
+      }
+    }
     
     if ('indiff' %in% colnames(data)) {
       # Plot empirical indifference points
diff --git a/R/internal_utils.R b/R/internal_utils.R
index 15c0aec..bba2c09 100644
--- a/R/internal_utils.R
+++ b/R/internal_utils.R
@@ -34,3 +34,23 @@ get_candidate_discount_functions <- function(arg) {
   return(unique(candidates))
 }
 
+get_pimm_func <- function(mod) {
+  # Get a function to compute the probability of selecting the immediate reward
+  
+  if (is(mod, 'td_bcnm')) {
+    frame <- do.call(get_prob_mod_frame, mod$config)
+  } else if (is(mod, 'td_ddm')) {
+    linpred_func <- do.call(get_linpred_func_ddm, mod$config)
+    frame <- function(data, par) {
+      drift <- linpred_func(data, par)
+      return(pimm_ddm(drift, par))
+    }
+  }
+  # Get a function where the parameter values are fixed
+  func <- function(data) {
+    frame(data, coef(mod))
+  }
+  
+  return(func)
+  
+}
\ No newline at end of file
diff --git a/README.Rmd b/README.Rmd
index c99fdbc..ff6c7ef 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -103,7 +103,7 @@ We can explore an even wider range of discount functions using nonlinear modelin
 
 ```{r}
 mod <- td_bcnm(td_bc_single_ptpt, discount_function = 'all')
-plot(mod, log = 'x', verbose = F)
+plot(mod, log = 'x', verbose = F, p_lines = c(0.05, 0.95))
 ```
 
 ### Drift diffusion models
diff --git a/README.md b/README.md
index ed00740..b2ad952 100644
--- a/README.md
+++ b/README.md
@@ -215,7 +215,7 @@ all of the following models are tested and the best-fitting one
 
 ``` r
 mod <- td_bcnm(td_bc_single_ptpt, discount_function = 'all')
-plot(mod, log = 'x', verbose = F)
+plot(mod, log = 'x', verbose = F, p_lines = c(0.05, 0.95))
 ```
 
 <img src="man/figures/README-unnamed-chunk-13-1.png" width="100%" />
diff --git a/man/figures/README-unnamed-chunk-13-1.png b/man/figures/README-unnamed-chunk-13-1.png
index 38389d9..c4dc992 100644
Binary files a/man/figures/README-unnamed-chunk-13-1.png and b/man/figures/README-unnamed-chunk-13-1.png differ
diff --git a/man/most_consistent_indiffs.Rd b/man/most_consistent_indiffs.Rd
new file mode 100644
index 0000000..2c25abc
--- /dev/null
+++ b/man/most_consistent_indiffs.Rd
@@ -0,0 +1,14 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{most_consistent_indiffs}
+\alias{most_consistent_indiffs}
+\title{Experimental method for computing indifference points}
+\usage{
+most_consistent_indiffs(data)
+}
+\arguments{
+\item{data}{Responses to score}
+}
+\description{
+Experimental method for computing indifference points
+}
diff --git a/man/plot.td_um.Rd b/man/plot.td_um.Rd
index 0ea5726..280d690 100644
--- a/man/plot.td_um.Rd
+++ b/man/plot.td_um.Rd
@@ -8,6 +8,8 @@
   x,
   type = c("summary", "endpoints", "link", "rt"),
   legend = T,
+  p_lines = NULL,
+  p_tol = 0.001,
   verbose = T,
   del = NULL,
   val_del = NULL,
@@ -22,6 +24,10 @@
 
 \item{legend}{Logical: display a legend? Only relevant for \code{type = 'summary'} and \code{type = 'rt'}}
 
+\item{p_lines}{Numerical vector. When \code{type = 'summary'} the discount curve, where the probability of selecting the immediate reward is 0.5, is plotted. \code{p_lines} allows you to specify other probabilities for which similar curves should be plotted (only applicable for probabilistic models, e.g. \code{td_bcnm}, \code{td_bclm} and \code{td_ddm})}
+
+\item{p_tol}{If \code{p_lines} is not \code{NULL}, what is the maximum distance that estimated probabilities can be from their true values? Smaller values results in slower plot generation}
+
 \item{verbose}{Whether to print info about, e.g., setting del = ED50 when \code{type = 'endpoints'}}
 
 \item{del}{Plots data for a particular delay}
@@ -39,7 +45,7 @@ Plot delay discounting models
 \dontrun{
 data("td_bc_single_ptpt")
 mod <- td_bclm(td_bc_single_ptpt, model = 'hyperbolic.1')
-plot(mod, type = 'summary')
+plot(mod, type = 'summary', p_lines = c(0.25, 0.75), log = 'x')
 plot(mod, type = 'endpoints')
 }
 }
diff --git a/tests/testthat/test-td_bclm.R b/tests/testthat/test-td_bclm.R
index 7a7e17e..f47b5be 100644
--- a/tests/testthat/test-td_bclm.R
+++ b/tests/testthat/test-td_bclm.R
@@ -34,6 +34,7 @@ while (model_idx <= length(models)) {
   test_that('plots', {
     expect_no_error(plot(mod, type = 'summary', verbose = F))
     expect_no_error(plot(mod, type = 'summary', verbose = F, log = 'x'))
+    expect_no_error(plot(mod, type = 'summary', verbose = F, log = 'x', p_lines = 0.1, p_tol = 0.1))
     expect_no_error(plot(mod, type = 'endpoints', verbose = F))
     expect_output(plot(mod, type = 'summary', verbose = T))
     expect_output(plot(mod, type = 'endpoints', verbose = T))
diff --git a/tests/testthat/test-td_bcnm.R b/tests/testthat/test-td_bcnm.R
index dff235f..4462af1 100644
--- a/tests/testthat/test-td_bcnm.R
+++ b/tests/testthat/test-td_bcnm.R
@@ -68,6 +68,7 @@ while (arg_combo_idx <= nrow(arg_combos)) {
   test_that('plots', {
     expect_no_error(plot(mod, type = 'summary', verbose = F))
     expect_no_error(plot(mod, type = 'summary', verbose = F, log = 'x'))
+    expect_no_error(plot(mod, type = 'summary', verbose = F, log = 'x', p_lines = 0.1, p_tol = 0.1))
     expect_no_error(plot(mod, type = 'endpoints', verbose = F))
     expect_output(plot(mod, type = 'endpoints', verbose = T))
     expect_no_error(plot(mod, type = 'endpoints', verbose = F, del = 100, val_del = 50))
diff --git a/tests/testthat/test-td_ddm.R b/tests/testthat/test-td_ddm.R
index 15477d7..32e4cbb 100644
--- a/tests/testthat/test-td_ddm.R
+++ b/tests/testthat/test-td_ddm.R
@@ -41,6 +41,7 @@ test_that('plots', {
   pdf(NULL) # Don't actually produce plots
   expect_no_error(plot(mod, type = 'summary', verbose = F))
   expect_no_error(plot(mod, type = 'summary', verbose = F, log = 'x'))
+  expect_no_error(plot(mod, type = 'summary', verbose = F, log = 'x', p_lines = 0.1, p_tol = 0.1))
   expect_no_error(plot(mod, type = 'endpoints', verbose = F))
   expect_output(plot(mod, type = 'summary', verbose = T))
   expect_output(plot(mod, type = 'endpoints', verbose = T))
diff --git a/tests/testthat/test-td_ipm.R b/tests/testthat/test-td_ipm.R
index 887ea00..9a40d53 100644
--- a/tests/testthat/test-td_ipm.R
+++ b/tests/testthat/test-td_ipm.R
@@ -80,4 +80,5 @@ test_that('errors', {
   expect_error(td_ipm(data.frame(del = 1:10)))
   expect_error(td_ipm(df, discount_function = 'new'))
   expect_error(plot(mod, type = 'endpoints'))
+  expect_error(plot(mod, p_lines = 0.1))
 })
\ No newline at end of file
diff --git a/vignettes/visualizing-models.Rmd b/vignettes/visualizing-models.Rmd
index 015a321..5a9b1b1 100644
--- a/vignettes/visualizing-models.Rmd
+++ b/vignettes/visualizing-models.Rmd
@@ -36,6 +36,12 @@ The plotting function prints some info, telling us it is plotting the discount c
 plot(mod, type = 'summary', verbose = F, log = 'x')
 ```
 
+Additionally, we can plot some information about how stochastic the individual's decision making was. The discount curve shows where the probability of selecting the immediate reward is predicted to be 50%, but we can plot curves for other probabilities as well. For example, we can show where the probability of selecting the immediate reward is 10% and 90% by setting `p_lines = c(0.1, 0.9)`. For more stochastic decision makers, there will be a greater separation between these (note that this is *not* a confidence interval for the discount curve itself):
+
+```{r}
+plot(mod, type = 'summary', verbose = F, log = 'x', p_lines = c(0.1, 0.9))
+```
+
 For an indifference point model, the discount function is usually plotted alongside the empirical indifference points:
 
 ```{r}