Here you will develop "d,p,q,r" functions for a certain distribution family. We'll call the family "accum" for "accumulate." The setting is that of repeatedly rolling a pair of dice. The random variable X is the number of rolls needed to achieve an accumulated total of at least k dots. So for instance the support of X ranges from ceiling(k/12) t…

accum.R

# P(dots == n) for k dice
dkdice <- function(n, k=2) {
  if (k == 1) {
    if (n >= 1 && n <= 6) {
      return(1/6)
    } else {
      return(0);
    }
  } else if (k < 1) {
    return(0);
  }
  P <- 0
  for (d in 1:6) {
    P <- P + 1/6 * dkdice(n-d, k-1)
  }
  P
}

# P(dots <= n) for k dice
pkdice <- function(n, k=2) {
  P <- 0
  for (i in 1:n) {
    P <- P + dkdice(i, k)
  }
  return(P);
}

# Probability for needing i rolls to get at least k dots.
daccum <- function(i, k) {
  if (k > 0 && i == 1) {
    return(1-pkdice(k-1, 2));
  } else if (i == 0) {
    return(0);
  }
  P <- 0
  for (roll in 1:12) {
    if (roll > k) {
      break;
    }
    # Assume we roll a `roll` in this round; chain with P(dots >= k-roll on i-1 remaining rolls)
    P <- P + daccum(i-1, k-roll) * dkdice(roll, 2)
  }
  P
}

# Check that probabilities sum to 1 for all k's up to a certain limit.
test_daccum <- function(upto_k) {
  for (k in 2:upto_k) {
    mini <- ceiling(k/12)
    maxi <- ceiling(k/2)

    p <- 0
    for (i in mini:maxi) {
      p <- p + daccum(i, k)
    }
    cat(p, "\n")
  }
}

# Probability for needing i or fewer rolls to get at least k dots.
paccum <- function(i, k) {
  if (i < 1) {
    return(0);
  }
  P <- 0
  for (rolls in 1:i) {
    P <- P + daccum(rolls, k)
  }
  P
}

# Find number of rolls i such that P(X = i for >= k dots) >= q
qaccum <- function(q, k) {
  i <- 1
  while (daccum(i,k) < q) {
    i <- i+1
  }
  i
}

raccum <- function(nreps, k) {
  s <- rep(0, nreps)
  pvector <- mapply(dkdice, 1:12)
  dots <- 1:12
  for (i in 1:nreps) {
    n <- 0
    j <- 0
    # Roll dice until accumulated number of dots equals or exceeds k.
    while (n < k) {
      n <- n + sample(dots, 1, prob=pvector)
      j <- j+1
    }
    s[i] <- j
  }
  s
}