diff --git a/DESCRIPTION b/DESCRIPTION
index c9141c06..9a2c9657 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -40,3 +40,6 @@ Imports:
 URL: https://github.com/venpopov/bmm, https://venpopov.github.io/bmm/
 BugReports: https://github.com/venpopov/bmm/issues
 VignetteBuilder: knitr
+Depends: 
+    R (>= 2.10)
+LazyData: true
diff --git a/NAMESPACE b/NAMESPACE
index fdeaddb2..35d53339 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -24,6 +24,9 @@ export(c_sqrtexp2bessel)
 export(calc_error_relative_to_nontargets)
 export(check_data)
 export(configure_model)
+export(dIMM)
+export(dmixture2p)
+export(dmixture3p)
 export(dsdm)
 export(fit_model)
 export(gen_3p_data)
@@ -32,10 +35,19 @@ export(get_model_prior)
 export(k2sd)
 export(mixture2p)
 export(mixture3p)
+export(pIMM)
+export(pmixture2p)
+export(pmixture3p)
 export(postprocess_brm)
 export(print_pretty_models_md)
 export(psdm)
+export(qIMM)
+export(qmixture2p)
+export(qmixture3p)
 export(qsdm)
+export(rIMM)
+export(rmixture2p)
+export(rmixture3p)
 export(rsdm)
 export(sdmSimple)
 export(softmax)
diff --git a/R/data.R b/R/data.R
new file mode 100644
index 00000000..fa642885
--- /dev/null
+++ b/R/data.R
@@ -0,0 +1,45 @@
+#' Data from Experiment 2 reported by Zhang & Luck (2008)
+#'
+#' Raw data of 8 subjects for the response error in a continuous reproduction task
+#' with set size 1, 2, 3, and 6 reported by Zhang & Luck (2008).
+#'
+#' @format ## `ZhangLuck_2008`
+#' A data frame with 4,000 rows and 9 columns:
+#' \describe{
+#'   \item{subID}{Integer uniquely identifying different subjects}
+#'   \item{trial}{Trial identifyier}
+#'   \item{setsize}{The setsize of the data in this row}
+#'   \item{RespErr}{The response error, that is the difference between the reponse
+#'   given and the target color.}
+#'   \item{Pos_Lure1, Pos_Lure2, Pos_Lure3, Pos_Lure4, Pos_Lure5}{Position of the lure items relative to the target color.}
+#'
+#' }
+#'
+#' @source <https://www.nature.com/articles/nature06860>
+"ZhangLuck_2008"
+
+
+#' Data from Experiment 1 reported by Oberauer & Lin (2017)
+#'
+#' Raw data of 19 subjects that completed a continuous reproduction task
+#' with set size 1 to 8 reported by Oberauer & Lin (2017).
+#'
+#' @format ## `OberauerLin_2017`
+#' A data frame with 4,000 rows and 9 columns:
+#' \describe{
+#'   \item{ID}{Integer uniquely identifying different subjects}
+#'   \item{Session}{Session number}
+#'   \item{Trial}{Trial number within each session}
+#'   \item{SetSize}{The setsize of the data in this row}
+#'   \item{Response}{The response in degrees given on the color wheel}
+#'   \item{deviation}{The response error or deviation of the `Response` from the target color (i.e., `Item1_Col`) in degrees.}
+#'   \item{dev_rad}{The response error converted from degrees to radians.}
+#'   \item{Item1_Col,Item2_Col,Item3_Col,Item4_Col,Item5_Col,Item6_Col,Item7_Col,Item8_Col}{The absolute colors of all items in degrees of the color wheel. Although there are always eight values given even for set sizes smaller than 8, only the colors from item 1 to the respective set size were shown.}
+#'   \item{Item1_Pos,Item2_Pos,Item3_Pos,Item4_Pos,Item5_Pos,Item6_Pos,Item7_Pos,Item8_Pos}{The position of all items in clockwise order around the circle. There were 12 possible positions, thus each position was 30 degrees apart from each other. Although positions are always given for all items, only item 1 to the respective set size was shown.}
+#'   \item{Item1_Col_rad,Item2_Col_rad,Item3_Col_rad,Item4_Col_rad,Item5_Col_rad,Item6_Col,Item7_Col_rad,Item8_Col_rad}{The relative position of colors to the target item (i.e. `Item1_Col`) of all items in radians. Although there are always eight values given even for set sizes smaller than 8, only the colors from item 1 to the respective set size were shown.}
+#'   \item{Item1_Pos_rad,Item2_Pos_rad,Item3_Pos_rad,Item4_Pos_rad,Item5_Pos_rad,Item6_Pos,Item7_Pos_rad,Item8_Pos_rad}{The relative position of all items to the position of the target item (i.e. `Item1_Pos`) in radians. Although there are always eight values given even for set sizes smaller than 8, only the colors from item 1 to the respective set size were shown.}
+#'
+#' }
+#'
+#' @source <https://osf.io/m4shu>
+"OberauerLin_2017"
diff --git a/R/distributions.R b/R/distributions.R
index 6c84f14d..b3b63416 100644
--- a/R/distributions.R
+++ b/R/distributions.R
@@ -160,7 +160,6 @@ qsdm <- function(p, mu=0, c=3, kappa=3.5, parametrization = "sqrtexp") {
   .NotYetImplemented()
 }
 
-
 #' @rdname SDMdist
 #' @export
 rsdm <- function(n, mu = 0, c = 3, kappa = 3.5, parametrization = "sqrtexp") {
@@ -195,9 +194,6 @@ rsdm <- function(n, mu = 0, c = 3, kappa = 3.5, parametrization = "sqrtexp") {
   .rsdm_inner(n, mu, c, kappa, parametrization, xa)
 }
 
-
-
-
 # helper functions for calculating the density of the SDM distribution
 .dsdm_numer_bessel <- function(x, mu, c, kappa, log=FALSE) {
   if (isTRUE(any(kappa < 0))) {
@@ -241,4 +237,368 @@ rsdm <- function(n, mu = 0, c = 3, kappa = 3.5, parametrization = "sqrtexp") {
 .dsdm_integrate_numer_v <- Vectorize(.dsdm_integrate_numer,
                                      vectorize.args = c("mu", "c", "kappa", 'lower','upper'))
 
+#' @title The two-parameter mixture model (mixture2p)
+#'
+#' @description Density, distribution, and random generation functions for the
+#'   two-parameter mixture model with the location of `mu`, precision of memory
+#'   representations `kappa` and probability of recalling items from memory `pMem`.
+#'
+#' @name mixture2p_dist
+#'
+#' @param x Vector of observed responses
+#' @param q Vector of quantiles
+#' @param p Vector of probability
+#' @param n Number of observations to generate data for
+#' @param mu Vector of locations
+#' @param kappa Vector of precision values
+#' @param pMem Vector of probabilities for memory recall
+#' @param log Logical; if `TRUE`, values are returned on the log scale.
+#'
+#' @keywords distribution
+#'
+#' @references Zhang, W., & Luck, S. J. (2008). Discrete fixed-resolution
+#'   representations in visual working memory. Nature, 453.
+#'
+#' @return `dmixture2p` gives the density of the two-parameter mixture model, `pmixture2p`
+#'   gives the cumulative distribution function of the two-parameter mixture model,
+#'   `qmixture2p` gives the quantile function of the two-parameter mixture model, and
+#'   `rmixture2p` gives the random generation function for the two-parameter mixture model.
+#'
+#' @export
+#'
+#' @examples
+#' # example code
+#'
+dmixture2p <- function(x, mu=0, kappa=5, pMem = 0.6, log = FALSE) {
+  if (isTRUE(any(kappa < 0))) {
+    stop("kappa must be non-negative")
+  }
+
+  if (isTRUE(any(pMem < 0))) {
+    stop("pMem must be larger than zero.")
+  }
+
+  if (isTRUE(any(pMem > 1))) {
+    stop("pMem must be smaller than one.")
+  }
+
+  density <- matrix(data = NaN, nrow = length(x), ncol = 2)
+
+  density[,1] <- log(pMem) + brms::dvon_mises(x = x,mu = mu , kappa = kappa, log = T)
+  density[,2] <- log(1 - pMem) + brms::dvon_mises(x = x,mu = 0 , kappa = 0, log = T)
+
+  density <- matrixStats::rowLogSumExps(density)
+
+  if (log == FALSE) {
+    return(exp(density))
+  }
+
+  return(density)
+}
+
+#' @rdname mixture2p_dist
+#' @export
+pmixture2p <- function(q, mu=0, kappa=7, pMem = 0.8) {
+  .NotYetImplemented()
+}
+
+#' @rdname mixture2p_dist
+#' @export
+qmixture2p <- function(p, mu=0, kappa=5, pMem = 0.6) {
+  .NotYetImplemented()
+}
+
+#' @rdname mixture2p_dist
+#' @export
+rmixture2p <- function(n, mu=0, kappa=5, pMem = 0.6) {
+  if (isTRUE(any(kappa < 0))) {
+    stop("kappa must be non-negative")
+  }
+
+  if (isTRUE(any(pMem < 0))) {
+    stop("pMem must be larger than zero.")
+  }
+
+  if (isTRUE(any(pMem > 1))) {
+    stop("pMem must be smaller than one.")
+  }
+
+  maxy <- dmixture2p(0, 0, kappa, pMem)
+  xa <- c()
+
+  .rmixture2p_inner <- function(n, mu, c, kappa, pMem, xa) {
+    x <- stats::runif(n, -pi, pi)
+    y <- stats::runif(n, 0, 1) * maxy
+    accept <- y < dmixture2p(x, mu, kappa, pMem)
+    xa <- c(xa, x[accept])
+
+    if (length(xa) < n) {
+      return(.rmixture2p_inner(n, mu, c, kappa, pMem, xa))
+    }
+
+    return(xa[1:n])
+  }
+
+  .rmixture2p_inner(n, mu, c, kappa, pMem, xa)
+}
+
+
+#' @title The three-parameter mixture model (mixture3p)
+#'
+#' @description Density, distribution, and random generation functions for the
+#'   three-parameter mixture model with the location of `mu`, precision of memory
+#'   representations `kappa`, probability of recalling items from memory `pMem`,
+#'   and probability of recalling non-targets `pNT`.
+#'
+#' @name mixture3p_dist
+#'
+#' @param x Vector of observed responses
+#' @param q Vector of quantiles
+#' @param p Vector of probability
+#' @param n Number of observations to generate data for
+#' @param mu Vector of locations. First value represents the location of the target item
+#'   and any additional values indicate the location of non-target items.
+#' @param kappa Vector of precision values
+#' @param pMem Vector of probabilities for memory recall
+#' @param pNT Vector of probabilities for swap errors
+#' @param log Logical; if `TRUE`, values are returned on the log scale.
+#'
+#' @keywords distribution
+#'
+#' @references Bays, P. M., Catalao, R. F. G., & Husain, M. (2009). The precision of visual
+#'   working memory is set by allocation of a shared resource. Journal of Vision, 9(10), 7.
+#'
+#' @return `dmixture3p` gives the density of the three-parameter mixture model, `pmixture3p`
+#'   gives the cumulative distribution function of the two-parameter mixture model,
+#'   `qmixture3p` gives the quantile function of the two-parameter mixture model, and
+#'   `rmixture3p` gives the random generation function for the two-parameter mixture model.
+#'
+#' @export
+#'
+#' @examples
+#' # example code
+#'
+dmixture3p <- function(x, mu=c(0,2,-1.5), kappa = 5, pMem = 0.6, pNT = 0.2, log = FALSE) {
+  if (isTRUE(any(kappa < 0))) {
+    stop("kappa must be non-negative")
+  }
+
+  if (isTRUE(any(pMem < 0))) {
+    stop("pMem must be larger than zero.")
+  }
+
+  if (isTRUE(any(pNT < 0))) {
+    stop("pMem must be larger than zero.")
+  }
 
+  if (isTRUE(any(pMem + pNT > 1))) {
+    stop("The sum of pMem and pNT must be smaller than one.")
+  }
+
+  density <- matrix(data = NaN, nrow = length(x), ncol = length(mu) + 1)
+  probs <- c(pMem,
+             rep(pNT/(length(mu) - 1), each = length(mu) - 1),
+             (1 - pMem - pNT))
+
+  for (i in 1:(length(mu))) {
+    density[,i] <- log(probs[i]) +
+      brms::dvon_mises(x = x, mu = mu[i], kappa = kappa, log = T)
+  }
+
+  density[,length(mu) + 1] <- log(probs[length(mu) + 1]) +
+    stats::dunif(x = x,-pi, pi, log = T)
+
+  density <- matrixStats::rowLogSumExps(density)
+
+  if (log == FALSE) {
+    return(exp(density))
+  }
+
+  # return the weighted density
+  return(density)
+}
+
+#' @rdname mixture3p_dist
+#' @export
+pmixture3p <- function(q, mu=c(0,2,-1.5), kappa = 5, pMem = 0.6, pNT = 0.2) {
+  .NotYetImplemented()
+}
+
+#' @rdname mixture3p_dist
+#' @export
+qmixture3p <- function(p, mu=c(0,2,-1.5), kappa = 5, pMem = 0.6, pNT = 0.2) {
+  .NotYetImplemented()
+}
+
+#' @rdname mixture3p_dist
+#' @export
+rmixture3p <- function(n, mu=c(0,2,-1.5), kappa = 5, pMem = 0.6, pNT = 0.2) {
+  if (isTRUE(any(kappa < 0))) {
+    stop("kappa must be non-negative")
+  }
+
+  if (isTRUE(any(pMem < 0))) {
+    stop("pMem must be larger than zero.")
+  }
+
+  if (isTRUE(any(pNT < 0))) {
+    stop("pMem must be larger than zero.")
+  }
+
+  if (isTRUE(any(pMem + pNT > 1))) {
+    stop("The sum of pMem and pNT must be smaller than one.")
+  }
+
+  maxy <- dmixture3p(mu[1], mu, kappa, pMem, pNT)
+  xa <- c()
+
+  .rmixture3p_inner <- function(n, mu, c, kappa, pMem, pNT, xa) {
+    x <- stats::runif(n, -pi, pi)
+    y <- stats::runif(n, 0, 1) * maxy
+    accept <- y < dmixture3p(x, mu, kappa, pMem, pNT)
+    xa <- c(xa, x[accept])
+
+    if (length(xa) < n) {
+      return(.rmixture3p_inner(n, mu, c, kappa, pMem, pNT, xa))
+    }
+
+    return(xa[1:n])
+  }
+
+  .rmixture3p_inner(n, mu, c, kappa, pMem, pNT, xa)
+}
+
+#' @title The Interference Measurement Model (IMM)
+#'
+#' @description Density, distribution, and random generation functions for the
+#'   interference measurement model with the location of `mu`, strength of cue-
+#'   dependent activation `c`, strength of cue-independent activation `a`, the
+#'   generalization gradient `s`, and the precision of memory representations
+#'   `kappa`.
+#'
+#' @name IMMdist
+#'
+#' @param x Vector of observed responses
+#' @param q Vector of quantiles
+#' @param p Vector of probability
+#' @param n Number of observations to generate data for
+#' @param mu Vector of locations
+#' @param dist Vector of distances of the item locations to the cued location
+#' @param kappa Vector of precision values
+#' @param c Vector of strengths for cue-dependent activation
+#' @param a Vector of strengths for cue-independent activation
+#' @param s Vector of generalization gradients
+#' @param b Vector of baseline activation
+#' @param log Logical; if `TRUE`, values are returned on the log scale.
+#'
+#' @keywords distribution
+#'
+#' @references Oberauer, K., Stoneking, C., Wabersich, D., & Lin, H.-Y. (2017).
+#'   Hierarchical Bayesian measurement models for continuous reproduction of visual
+#'   features from working memory. Journal of Vision, 17(5), 11.
+#'
+#' @return `dIMM` gives the density of the interference measurement model, `pIMM`
+#'   gives the cumulative distribution function of the interference measurement model,
+#'   `qIMM` gives the quantile function of the interference measurement model, and
+#'   `rIMM` gives the random generation function for the interference measurement model.
+#'
+#' @export
+#'
+#' @examples
+#' # example code
+#'
+dIMM <- function(x, mu=c(0,2,-1.5), dist = c(0,0.5,2),
+                 c=1, a = 0.2, b = 0, s = 2, kappa=5, log = FALSE) {
+  if (isTRUE(any(kappa < 0))) {
+    stop("kappa must be non-negative")
+  }
+
+  if (length(mu) != length(dist)) {
+    stop("The number of items does not match the distances provided from the cued location.")
+  }
+
+  if (isTRUE(any(s < 0))) {
+    stop("s must be non-negative")
+  }
+
+  if (isTRUE(any(dist < 0))) {
+    stop("all distances have to be positive.")
+  }
+
+  # compute activation for all items
+  acts <- rep(c, length(mu)) * exp(-s*dist) + rep(a, length(mu))
+
+  # add activation of background noise
+  acts <- c(acts,b)
+
+  # compute probability for responding stemming from each distribution
+  probs <- exp(acts)/sum(exp(acts))
+
+  density <- matrix(data = NaN, nrow = length(x), ncol = length(mu) + 1)
+
+  for (i in 1:(length(mu))) {
+    density[,i] <- log(probs[i]) +
+      brms::dvon_mises(x = x, mu = mu[i], kappa = kappa, log = T)
+  }
+
+  density[,length(mu) + 1] <- log(probs[length(mu) + 1]) +
+    stats::dunif(x = x,-pi, pi, log = T)
+
+  density <- matrixStats::rowLogSumExps(density)
+
+  if (log == FALSE) {
+    return(exp(density))
+  }
+
+  # return the weighted density
+  return(density)
+}
+
+#' @rdname IMMdist
+#' @export
+pIMM <- function(q, mu=c(0,2,-1.5), dist = c(0,0.5,2),
+                 c=1, a = 0.2, b = 0, s = 2, kappa=5) {
+  .NotYetImplemented()
+}
+
+#' @rdname IMMdist
+#' @export
+qIMM <- function(p, mu=c(0,2,-1.5), dist = c(0,0.5,2),
+                 c=1, a = 0.2, b = 0, s = 2, kappa=5) {
+  .NotYetImplemented()
+}
+
+#' @rdname IMMdist
+#' @export
+rIMM <- function(n, mu=c(0,2,-1.5), dist = c(0,0.5,2),
+                 c=1, a = 0.2, b = 0, s = 2, kappa=5) {
+  if (isTRUE(any(kappa < 0))) {
+    stop("kappa must be non-negative")
+  }
+
+  if (length(mu) != length(dist)) {
+    stop("The number of items does not match the distances provided from the cued location.")
+  }
+
+  if (isTRUE(any(s < 0))) {
+    stop("s must be non-negative")
+  }
+
+  maxy <- dIMM(mu[1], mu, dist, c, a, b, s, kappa)
+  xa <- c()
+
+  .rIMM_inner <- function(n, mu, dist, c, a, b, s, kappa, xa) {
+    x <- stats::runif(n, -pi, pi)
+    y <- stats::runif(n, 0, 1) * maxy
+    accept <- y < dIMM(x, mu, dist, c, a, b, s, kappa)
+    xa <- c(xa, x[accept])
+
+    if (length(xa) < n) {
+      return(.rIMM_inner(n, mu, dist, c, a, b, s, kappa, xa))
+    }
+
+    return(xa[1:n])
+  }
+
+  .rIMM_inner(n, mu, dist, c ,a ,b ,s ,kappa , xa)
+}
diff --git a/data/OberauerLin_2017.rda b/data/OberauerLin_2017.rda
new file mode 100644
index 00000000..3b3e4c13
Binary files /dev/null and b/data/OberauerLin_2017.rda differ
diff --git a/data/ZhangLuck_2008.rda b/data/ZhangLuck_2008.rda
new file mode 100644
index 00000000..ae142ad5
Binary files /dev/null and b/data/ZhangLuck_2008.rda differ
diff --git a/inst/imm_vignette_fit.rds b/inst/imm_vignette_fit.rds
new file mode 100644
index 00000000..67f25843
Binary files /dev/null and b/inst/imm_vignette_fit.rds differ
diff --git a/man/IMMdist.Rd b/man/IMMdist.Rd
new file mode 100644
index 00000000..54400aff
--- /dev/null
+++ b/man/IMMdist.Rd
@@ -0,0 +1,103 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/distributions.R
+\name{IMMdist}
+\alias{IMMdist}
+\alias{dIMM}
+\alias{pIMM}
+\alias{qIMM}
+\alias{rIMM}
+\title{The Interference Measurement Model (IMM)}
+\usage{
+dIMM(
+  x,
+  mu = c(0, 2, -1.5),
+  dist = c(0, 0.5, 2),
+  c = 1,
+  a = 0.2,
+  b = 0,
+  s = 2,
+  kappa = 5,
+  log = FALSE
+)
+
+pIMM(
+  q,
+  mu = c(0, 2, -1.5),
+  dist = c(0, 0.5, 2),
+  c = 1,
+  a = 0.2,
+  b = 0,
+  s = 2,
+  kappa = 5
+)
+
+qIMM(
+  p,
+  mu = c(0, 2, -1.5),
+  dist = c(0, 0.5, 2),
+  c = 1,
+  a = 0.2,
+  b = 0,
+  s = 2,
+  kappa = 5
+)
+
+rIMM(
+  n,
+  mu = c(0, 2, -1.5),
+  dist = c(0, 0.5, 2),
+  c = 1,
+  a = 0.2,
+  b = 0,
+  s = 2,
+  kappa = 5
+)
+}
+\arguments{
+\item{x}{Vector of observed responses}
+
+\item{mu}{Vector of locations}
+
+\item{dist}{Vector of distances of the item locations to the cued location}
+
+\item{c}{Vector of strengths for cue-dependent activation}
+
+\item{a}{Vector of strengths for cue-independent activation}
+
+\item{b}{Vector of baseline activation}
+
+\item{s}{Vector of generalization gradients}
+
+\item{kappa}{Vector of precision values}
+
+\item{log}{Logical; if \code{TRUE}, values are returned on the log scale.}
+
+\item{q}{Vector of quantiles}
+
+\item{p}{Vector of probability}
+
+\item{n}{Number of observations to generate data for}
+}
+\value{
+\code{dIMM} gives the density of the interference measurement model, \code{pIMM}
+gives the cumulative distribution function of the interference measurement model,
+\code{qIMM} gives the quantile function of the interference measurement model, and
+\code{rIMM} gives the random generation function for the interference measurement model.
+}
+\description{
+Density, distribution, and random generation functions for the
+interference measurement model with the location of \code{mu}, strength of cue-
+dependent activation \code{c}, strength of cue-independent activation \code{a}, the
+generalization gradient \code{s}, and the precision of memory representations
+\code{kappa}.
+}
+\examples{
+# example code
+
+}
+\references{
+Oberauer, K., Stoneking, C., Wabersich, D., & Lin, H.-Y. (2017).
+Hierarchical Bayesian measurement models for continuous reproduction of visual
+features from working memory. Journal of Vision, 17(5), 11.
+}
+\keyword{distribution}
diff --git a/man/OberauerLin_2017.Rd b/man/OberauerLin_2017.Rd
new file mode 100644
index 00000000..7d5ac035
--- /dev/null
+++ b/man/OberauerLin_2017.Rd
@@ -0,0 +1,37 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data.R
+\docType{data}
+\name{OberauerLin_2017}
+\alias{OberauerLin_2017}
+\title{Data from Experiment 1 reported by Oberauer & Lin (2017)}
+\format{
+\subsection{\code{OberauerLin_2017}}{
+
+A data frame with 4,000 rows and 9 columns:
+\describe{
+\item{ID}{Integer uniquely identifying different subjects}
+\item{Session}{Session number}
+\item{Trial}{Trial number within each session}
+\item{SetSize}{The setsize of the data in this row}
+\item{Response}{The response in degrees given on the color wheel}
+\item{deviation}{The response error or deviation of the \code{Response} from the target color (i.e., \code{Item1_Col}) in degrees.}
+\item{dev_rad}{The response error converted from degrees to radians.}
+\item{Item1_Col,Item2_Col,Item3_Col,Item4_Col,Item5_Col,Item6_Col,Item7_Col,Item8_Col}{The absolute colors of all items in degrees of the color wheel. Although there are always eight values given even for set sizes smaller than 8, only the colors from item 1 to the respective set size were shown.}
+\item{Item1_Pos,Item2_Pos,Item3_Pos,Item4_Pos,Item5_Pos,Item6_Pos,Item7_Pos,Item8_Pos}{The position of all items in clockwise order around the circle. There were 12 possible positions, thus each position was 30 degrees apart from each other. Although positions are always given for all items, only item 1 to the respective set size was shown.}
+\item{Item1_Col_rad,Item2_Col_rad,Item3_Col_rad,Item4_Col_rad,Item5_Col_rad,Item6_Col,Item7_Col_rad,Item8_Col_rad}{The relative position of colors to the target item (i.e. \code{Item1_Col}) of all items in radians. Although there are always eight values given even for set sizes smaller than 8, only the colors from item 1 to the respective set size were shown.}
+\item{Item1_Pos_rad,Item2_Pos_rad,Item3_Pos_rad,Item4_Pos_rad,Item5_Pos_rad,Item6_Pos,Item7_Pos_rad,Item8_Pos_rad}{The relative position of all items to the position of the target item (i.e. \code{Item1_Pos}) in radians. Although there are always eight values given even for set sizes smaller than 8, only the colors from item 1 to the respective set size were shown.}
+
+}
+}
+}
+\source{
+\url{https://osf.io/m4shu}
+}
+\usage{
+OberauerLin_2017
+}
+\description{
+Raw data of 19 subjects that completed a continuous reproduction task
+with set size 1 to 8 reported by Oberauer & Lin (2017).
+}
+\keyword{datasets}
diff --git a/man/ZhangLuck_2008.Rd b/man/ZhangLuck_2008.Rd
new file mode 100644
index 00000000..80a6158c
--- /dev/null
+++ b/man/ZhangLuck_2008.Rd
@@ -0,0 +1,32 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data.R
+\docType{data}
+\name{ZhangLuck_2008}
+\alias{ZhangLuck_2008}
+\title{Data from Experiment 2 reported by Zhang & Luck (2008)}
+\format{
+\subsection{\code{ZhangLuck_2008}}{
+
+A data frame with 4,000 rows and 9 columns:
+\describe{
+\item{subID}{Integer uniquely identifying different subjects}
+\item{trial}{Trial identifyier}
+\item{setsize}{The setsize of the data in this row}
+\item{RespErr}{The response error, that is the difference between the reponse
+given and the target color.}
+\item{Pos_Lure1, Pos_Lure2, Pos_Lure3, Pos_Lure4, Pos_Lure5}{Position of the lure items relative to the target color.}
+
+}
+}
+}
+\source{
+\url{https://www.nature.com/articles/nature06860}
+}
+\usage{
+ZhangLuck_2008
+}
+\description{
+Raw data of 8 subjects for the response error in a continuous reproduction task
+with set size 1, 2, 3, and 6 reported by Zhang & Luck (2008).
+}
+\keyword{datasets}
diff --git a/man/mixture2p_dist.Rd b/man/mixture2p_dist.Rd
new file mode 100644
index 00000000..4781979e
--- /dev/null
+++ b/man/mixture2p_dist.Rd
@@ -0,0 +1,55 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/distributions.R
+\name{mixture2p_dist}
+\alias{mixture2p_dist}
+\alias{dmixture2p}
+\alias{pmixture2p}
+\alias{qmixture2p}
+\alias{rmixture2p}
+\title{The two-parameter mixture model (mixture2p)}
+\usage{
+dmixture2p(x, mu = 0, kappa = 5, pMem = 0.6, log = FALSE)
+
+pmixture2p(q, mu = 0, kappa = 7, pMem = 0.8)
+
+qmixture2p(p, mu = 0, kappa = 5, pMem = 0.6)
+
+rmixture2p(n, mu = 0, kappa = 5, pMem = 0.6)
+}
+\arguments{
+\item{x}{Vector of observed responses}
+
+\item{mu}{Vector of locations}
+
+\item{kappa}{Vector of precision values}
+
+\item{pMem}{Vector of probabilities for memory recall}
+
+\item{log}{Logical; if \code{TRUE}, values are returned on the log scale.}
+
+\item{q}{Vector of quantiles}
+
+\item{p}{Vector of probability}
+
+\item{n}{Number of observations to generate data for}
+}
+\value{
+\code{dmixture2p} gives the density of the two-parameter mixture model, \code{pmixture2p}
+gives the cumulative distribution function of the two-parameter mixture model,
+\code{qmixture2p} gives the quantile function of the two-parameter mixture model, and
+\code{rmixture2p} gives the random generation function for the two-parameter mixture model.
+}
+\description{
+Density, distribution, and random generation functions for the
+two-parameter mixture model with the location of \code{mu}, precision of memory
+representations \code{kappa} and probability of recalling items from memory \code{pMem}.
+}
+\examples{
+# example code
+
+}
+\references{
+Zhang, W., & Luck, S. J. (2008). Discrete fixed-resolution
+representations in visual working memory. Nature, 453.
+}
+\keyword{distribution}
diff --git a/man/mixture3p_dist.Rd b/man/mixture3p_dist.Rd
new file mode 100644
index 00000000..588a09fe
--- /dev/null
+++ b/man/mixture3p_dist.Rd
@@ -0,0 +1,66 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/distributions.R
+\name{mixture3p_dist}
+\alias{mixture3p_dist}
+\alias{dmixture3p}
+\alias{pmixture3p}
+\alias{qmixture3p}
+\alias{rmixture3p}
+\title{The three-parameter mixture model (mixture3p)}
+\usage{
+dmixture3p(
+  x,
+  mu = c(0, 2, -1.5),
+  kappa = 5,
+  pMem = 0.6,
+  pNT = 0.2,
+  log = FALSE
+)
+
+pmixture3p(q, mu = c(0, 2, -1.5), kappa = 5, pMem = 0.6, pNT = 0.2)
+
+qmixture3p(p, mu = c(0, 2, -1.5), kappa = 5, pMem = 0.6, pNT = 0.2)
+
+rmixture3p(n, mu = c(0, 2, -1.5), kappa = 5, pMem = 0.6, pNT = 0.2)
+}
+\arguments{
+\item{x}{Vector of observed responses}
+
+\item{mu}{Vector of locations. First value represents the location of the target item
+and any additional values indicate the location of non-target items.}
+
+\item{kappa}{Vector of precision values}
+
+\item{pMem}{Vector of probabilities for memory recall}
+
+\item{pNT}{Vector of probabilities for swap errors}
+
+\item{log}{Logical; if \code{TRUE}, values are returned on the log scale.}
+
+\item{q}{Vector of quantiles}
+
+\item{p}{Vector of probability}
+
+\item{n}{Number of observations to generate data for}
+}
+\value{
+\code{dmixture3p} gives the density of the three-parameter mixture model, \code{pmixture3p}
+gives the cumulative distribution function of the two-parameter mixture model,
+\code{qmixture3p} gives the quantile function of the two-parameter mixture model, and
+\code{rmixture3p} gives the random generation function for the two-parameter mixture model.
+}
+\description{
+Density, distribution, and random generation functions for the
+three-parameter mixture model with the location of \code{mu}, precision of memory
+representations \code{kappa}, probability of recalling items from memory \code{pMem},
+and probability of recalling non-targets \code{pNT}.
+}
+\examples{
+# example code
+
+}
+\references{
+Bays, P. M., Catalao, R. F. G., & Husain, M. (2009). The precision of visual
+working memory is set by allocation of a shared resource. Journal of Vision, 9(10), 7.
+}
+\keyword{distribution}
diff --git a/tests/testthat/test-distributions.R b/tests/testthat/test-distributions.R
index 86951902..7ac7c996 100644
--- a/tests/testthat/test-distributions.R
+++ b/tests/testthat/test-distributions.R
@@ -50,3 +50,56 @@ test_that("dsdm parametrization conversion returns accurate results", {
   d2 <- dsdm(y, 0, c_se, kappa, parametrization = "sqrtexp")
   expect_equal(d1,d2)
 })
+
+test_that("dmixture2p integrates to 1", {
+  expect_equal(integrate(dmixture2p, -pi,pi,
+                         mu = runif(1, min = -pi, pi),
+                         kappa = runif(1,min = 1, max = 20),
+                         pMem = runif(1, min = 0, max = 1))$value,1)
+})
+
+test_that("dmixture3p integrates to 1", {
+  expect_equal(integrate(dmixture3p, -pi,pi,
+                         mu = runif(3, min = -pi, pi),
+                         kappa = runif(1,min = 1, max = 20),
+                         pMem = runif(1, min = 0, max = 0.6),
+                         pNT = runif(1, min = 0, max = 0.3))$value,1)
+})
+
+test_that("dIMM integrates to 1", {
+  expect_equal(integrate(dIMM, -pi,pi,
+                         mu = runif(3, min = -pi, pi),
+                         dist = c(0, runif(2, min = 0.1, max = pi)),
+                         kappa = runif(1,min = 1, max = 20),
+                         c = runif(1, min = 0, max = 3),
+                         a = runif(1, min = 0, max = 1),
+                         s = runif(1, min = 1, max = 20),
+                         b = 0)$value,1)
+})
+
+test_that("rmixture2p returns values between -pi and pi", {
+  res <- rmixture2p(500, mu = runif(1, min = -pi, pi),
+              kappa = runif(1,min = 1, max = 20),
+              pMem = runif(1, min = 0, max = 1))
+  expect_true(all(res >= -pi) && all(res <= pi))
+})
+
+test_that("rmixture3p returns values between -pi and pi", {
+  res <- rmixture3p(500, mu = runif(3, min = -pi, pi),
+                    kappa = runif(1,min = 1, max = 20),
+                    pMem = runif(1, min = 0, max = 0.6),
+                    pNT = runif(1, min = 0, max = 0.3))
+  expect_true(all(res >= -pi) && all(res <= pi))
+})
+
+test_that("rIMM returns values between -pi and pi", {
+  res <- rIMM(500, mu = runif(3, min = -pi, pi),
+                    dist = c(0, runif(2, min = 0.1, max = pi)),
+                    kappa = runif(1,min = 1, max = 20),
+                    c = runif(1, min = 0, max = 3),
+                    a = runif(1, min = 0, max = 1),
+                    s = runif(1, min = 1, max = 20),
+                    b = 0)
+  expect_true(all(res >= -pi) && all(res <= pi))
+})
+
diff --git a/vignettes/IMM.Rmd b/vignettes/IMM.Rmd
new file mode 100644
index 00000000..25e23c97
--- /dev/null
+++ b/vignettes/IMM.Rmd
@@ -0,0 +1,250 @@
+---
+title: "The Interference Measurement Model (IMM)"
+output: bookdown::html_document2
+bibliography: REFERENCES.bib
+header-includes:
+  - \usepackage{amsmath}
+vignette: >
+  %\VignetteIndexEntry{The Interference Measurement Model (IMM)}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+pkgdown:
+  as_is: true
+---
+
+<style type="text/css">
+div.main-container {
+max-width: 800px !important;
+}
+
+p {
+margin-top: 1.5em ;
+margin-bottom: 1.5em ;
+}
+</style>
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+# Introduction to the model
+
+The Interference Measurement Model (IMM) is a measurement model for continous reproduction tasks in the domain of visual working memory. The model was introduced by @oberauerHierarchicalBayesianMeasurement2017 . The aim of the IMM is to capture the response behavior in continuous reproduction tasks including the occurrence of swap errors to other items encoded in visual working memory.
+
+The IMM assumes that during retrieval, to-be-remembered items or features (e.g., colors, orientations, or shapes) are associated with the context they appeared in (e.g. the spatial location). These associations can be continuous in their strength and represent bindings between the contents and context of to-be-remembered information (see Figure 1 of @oberauerHierarchicalBayesianMeasurement2017). At retrieval there are different sources of activation that contribute to the activation of the to-be-retrieved contents. Background noise (b) uniformly activates all possible responses, for example all 360 colors that participants can chose from in a color wheel experiment. Cue-independent activation (a) equally activates all features that were encoded into visual working memory during retrieval. And cue-dependent activation (c) activates the features that are associated with the current retrieval cue (e.g., the spatial location cued to be retrieved). Additionally, the IMM assumes that cue-dependent activation follows a generalization gradient (s) that activates similar contexts.
+
+The activation for each potential feature $x$ that could be retrieved is the sum of the weighted activation from all three activation sources, given a retrieval cue $L$ at the location $\theta$:
+
+$$
+A(x|L_\theta) = b \times A_b(x) + a \times A_a(x) + c \times A_c(c|L_\theta)
+$$
+
+The background activation ($A_b$) is independent from all encoded features and thus modeled as a uniform distribution around the circular feature space. This is implemented as a von-Mises (vM) distribution centered on 0\footnote{As this is a uniform distribution the location does not matter and could be fixed to an arbitrary value.} with a precision $\kappa = 0$:
+
+$$
+A_b(x) = vM(x,0,0)
+$$
+
+The cue-independent activation ($A_a$) is modeled as the sum of von Mises distributions centered on the feature values $x_i$ of all $n$ encoded items $i$ with the precision of memory \kappa:
+
+$$
+A_a(x) = \sum^n_{i = 1} vM(x,x_i,\kappa)
+$$
+
+The cue-dependent activation ($A_c$) is again modeled as the sum of von Mises distributions centered on the feature values $x_i$ of all $n$ encoded items $i$ with the precision of memory \kappa. These distributions are weighted by the spatial distance $D$ between the context $L$ a feature was associated with to the cue context $L_\theta$. This distance is weighted by the generalization gradient $s$ that captures the specificity of bindings to the cue dimension:
+
+$$
+A_c(x|L_\theta) = \sum^n_{i = 1} e^{-s*D(L,L\theta)} \times vM(x,x_i,\kappa)
+$$
+
+The probability for choosing each response $\hat{x}$ then results from normalizing the activation for all possible responses $N$. In the original publication of the IMM this was done using Luce's choice rule:
+
+$$
+P(\hat{x}|L_\theta) = \frac{A(\hat{x}|L_\theta)}{\sum^N_{j=1}A(j|L_\theta)}
+$$
+
+In the `bmm` package, we decided to implement an alternative normalization using the `softmax` normalization:
+
+$$
+P(\hat{x}|L_\theta) = \frac{e^{A(\hat{x}|L_\theta)}}{\sum^N_{j=1}e^{A(j|L_\theta)}}
+$$
+
+A comparison between these different normalization function in the context of activation based models of working memory can be found in the appendix of @oberauerSimpleMeasurementModels2019. Additionally, a more recent factorial comparison of different models for visual working memory @Oberauer_2023 indicated that the `softmax` normalization generally captures the observed data better than Luce's choice rule in the context of continuous reproduction tasks.
+
+In sum, the IMM assumes that responses in continuous reproduction tasks are the results of cue-based retrieval and cue-independent activation of all features corrupted by background noise.
+
+# Parametrization in the `bmm` package
+
+For identification, one of the weighting parameters has to be fixed. In the original publication the strenght of cue-dependent activation $c$ was fixed to one. The default setup of `brms` however currently only allows to fix the strength of the background noise $b$ to zero. Therefore, in all implementations of the IMM in the `bmm` package, the strength of cue-dependent and cue-independent activation, $c$ and $a$, can be estimated and predicted by independent variables.
+
+Apart from that, both the precision of memory representations $\kappa$ and the generalization gradient $s$ are parameterized the same way as in the original publication.
+
+Additionally, because we use the `softmax` normalization for translating activation into probabilities, the estimates for the strenght of cue-dependent and -independent activation, $c$ and $a$ have to be interpreted relatively to the strength of the baseline activation $b$ that is fixed to zero. Thus, it is possible that the strength of cue-dependent and cue-independent activation, $c$ and $a$, become negative. This does not reflect an absolute negative activation but rather an activation that is relatively smaller than the baseline activation. 
+
+# Fitting the model with `bmm`
+
+You should start by loading the `bmm` package:
+
+```{r}
+library(bmm)
+```
+
+## Generating simulated data
+
+Should you already have a data set you want to fit, you can skip this section. Alternatively, you can use data provided with the package (add reference to data) or generate data using the random generation function provided in the `bmm` package.
+
+```{r}
+# set seed for reproducibility
+set.seed(123)
+
+# specfiy generating parameters
+Cs <- c(4,4,2,2)
+As <- c(0.5,1,0.5,0.5)
+Ss <- c(10,10,5,5)
+kappas <- c(15,10,15,10)
+nTrials = 1000
+setsize = 5
+
+simData <- data.frame()
+for (i in 1:length(Cs)) {
+  # generate different non-target locations for each condition
+  item_location <- c(0, runif(setsize - 1, -pi,pi))
+  
+  # generate different distances for each condition
+  item_distance <- c(0, runif(setsize - 1, min = 0.1, max = pi))
+  
+  # simulate data for each condition
+  genData <- rIMM(n = nTrials,
+                  mu = item_location,
+                  dist = item_distance,
+                  c = Cs[i], a = As[i],
+                  b = 0, s = Ss[i], kappa = kappas[i])
+  
+  condData <- data.frame(
+    respErr = genData,
+    trialID = 1:nTrials,
+    cond = i,
+    color_item1 = 0,
+    dist_item1 = 0
+  )
+  
+  init_colnames <- colnames(condData)
+  
+  for (j in 1:(setsize - 1)) {
+    condData <- cbind(condData,item_location[j + 1])
+    condData <- cbind(condData,item_distance[j + 1])
+  }
+  
+  colnames(condData) <- c(init_colnames,
+                          paste0(rep(c("color_item","dist_item"),times = setsize - 1), 
+                                 rep(2:(setsize),each = 2)))
+  
+  simData <- rbind(simData,condData)
+}
+
+# convert condition variable to a factor
+
+simData$cond <- as.factor(simData$cond)
+```
+
+
+## Estimating the model with `bmm`
+
+To estimate the `IMM` we first need to specify a formula. In this formula, we directly estimate all parameters for each of the four conditions:
+
+```{r}
+model_formula <- brms::bf(respErr ~ 1,
+                          c ~ 0 + cond,
+                          a ~ 0 + cond,
+                          s ~ 0 + cond,
+                          kappa ~ 0 + cond)
+```
+
+Then, we can specify the model that we want to estimate:
+
+```{r}
+model <- IMMfull(non_target = paste0("color_item",2:setsize),
+                 setsize = setsize,
+                 spaPos = paste0("dist_item",2:setsize))
+```
+
+Finally, we can fit the model by passing all the relevant arguments to the `fit_model` function:
+
+``` r 
+fit <- fit_model(
+formula = model_formula,
+data = simData,
+model = model,
+parallel = TRUE,
+chains = 4,
+iter = 2000,
+backend = "cmdstanr"
+)
+```
+
+Running this model takes about 2 to 5 minutes (depending on the speed of your computer). Here we load an already saved fit object, but on your machine you will have to wait until the model finishes sampling. Both `brms` and `cmdstanr` typically print out some information on the sampling progress.
+
+```{r}
+fit <- readRDS(system.file("imm_vignette_fit.rds", package='bmm'))
+```
+
+Using this `fit` object we can have a quick look at the summary of the fitted model:
+
+```{r}
+summary(fit)
+```
+
+The first thing you might notice is that below the parts of the formula that was passed to the `fit_model` function, `bmm` has added a lot of additional specifications to implement the IMM. This is nothing that you have to check. But if you are interested in customizing and exploring different assumptions imposed by the IMM, you could start by taking this formula and adapting it accordingly.
+
+Next, we can have a look at the estimated parameters. The first thing we should check is if the sampling converged, this is indicated by all `Rhat` values being close to one. If you want to do more inspection of the sampling, you can check out the functionality implemented in `brms`to do this.
+
+The parameter estimates for `c` and `a` are already on their native scale, but both `s` and `kappa` are estimated using a `log` link function, so we have to transform these back to the native scale.
+
+```{r}
+fixedFX <- brms::fixef(fit)
+
+# print posterior means for the s parameter
+fixedFX[startsWith(rownames(fixedFX),"c_"),]
+
+# print posterior means for the s parameter
+fixedFX[startsWith(rownames(fixedFX),"a_"),]
+
+# print posterior means for the s parameter
+exp(fixedFX[grepl("s_",rownames(fixedFX)),])
+
+# print posterior means for the s parameter
+exp(fixedFX[grepl("kappa_",rownames(fixedFX)),])
+```
+
+These results indicate that all parameters, except for `s` were well recovered. As already noted by @oberauerHierarchicalBayesianMeasurement2017, a good recovery of the generalization gradient `s` requires a lot of data. Thus you might consider opting for the simplified version of the IMM without the `s` parameter, the `IMMabc`.
+
+We can illustrate the recovery of the data generating parameters by plotting the full posterior distributions alongside the data generating parameters. For this we need to extract the posterior draws using the `tidybayes` package and include the data generating parameters into the plots of the posteriors.
+
+```{r, message=FALSE, warning=FALSE}
+library(tidybayes)
+library(dplyr)
+library(tidyr)
+library(ggplot2)
+
+# extract the posterior draws
+draws <- tidybayes::tidy_draws(fit)
+draws <- select(draws, starts_with("b_")) %>% select(-(1:3)) %>% 
+  mutate_at(vars(starts_with("b_s")),exp) %>% 
+  mutate_at(vars(starts_with("b_kappa")),exp)
+
+# plot posterior with original parameters overlayed as diamonds
+as.data.frame(draws) %>% 
+  gather(par, value) %>% 
+  ggplot(aes(value, par)) +
+  tidybayes::stat_halfeyeh(normalize = "groups") +
+  geom_point(data = data.frame(par = colnames(draws),
+                               value = c(Cs,As,Ss,kappas)),
+             aes(value,par), color = "red",
+             shape = "diamond", size = 2.5) +
+  scale_x_continuous(lim=c(0,20))
+```
+
+# References
diff --git a/vignettes/REFERENCES.bib b/vignettes/REFERENCES.bib
index 1af09b19..817d8d33 100644
--- a/vignettes/REFERENCES.bib
+++ b/vignettes/REFERENCES.bib
@@ -71,3 +71,33 @@ @article{Oberauer_Lin_2017
 	doi          = {10.1037/rev0000044},
 	language     = {en}
 }
+
+@article{oberauerHierarchicalBayesianMeasurement2017,
+  title = {Hierarchical {{Bayesian}} Measurement Models for Continuous Reproduction of Visual Features from Working Memory},
+  author = {Oberauer, Klaus and Stoneking, Colin and Wabersich, Dominik and Lin, Hsuan-Yu},
+  year = {2017},
+  month = may,
+  journal = {Journal of Vision},
+  volume = {17},
+  number = {5},
+  pages = {11},
+  issn = {1534-7362},
+  doi = {10.1167/17.5.11},
+  urldate = {2022-10-10}
+}
+
+@article{oberauerSimpleMeasurementModels2019,
+  title = {Simple Measurement Models for Complex Working-Memory Tasks},
+  author = {Oberauer, Klaus and Lewandowsky, Stephan},
+  year = {2019},
+  month = nov,
+  journal = {Psychological Review},
+  volume = {126},
+  number = {6},
+  pages = {880--932},
+  publisher = {{American Psychological Association}},
+  issn = {0033-295X},
+  doi = {10.1037/rev0000159},
+  urldate = {2022-07-26}
+}
+
diff --git a/vignettes/sdm-simple.Rmd b/vignettes/sdm-simple.Rmd
index 62b2cdb2..37cb95d5 100644
--- a/vignettes/sdm-simple.Rmd
+++ b/vignettes/sdm-simple.Rmd
@@ -208,7 +208,7 @@ We can now inspect the results of the model fit:
 summary(fit)
 ```
 
-We see that all Rhat values are less than 1.01, which is a good sign that the chains have converged. In principle you have to do more inspection, but let us see the estimated parameters. The mmodel uses a log-link function for the `c` and `kappa` parameters, so we have to exponentiate the coefficients to get the estimated parameters:
+We see that all Rhat values are less than 1.01, which is a good sign that the chains have converged. In principle you have to do more inspection, but let us see the estimated parameters. The model uses a log-link function for the `c` and `kappa` parameters, so we have to exponentiate the coefficients to get the estimated parameters:
 
 ```{r}
 exp(brms::fixef(fit)[2:7,1])