diff --git a/DESCRIPTION b/DESCRIPTION
index f4fa1821..4ced1a8f 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: simtrial
 Type: Package
 Title: Clinical Trial Simulation
-Version: 0.3.0.9
+Version: 0.3.0.10
 Authors@R: c(
     person("Keaven", "Anderson", email = "keaven_anderson@merck.com", role = c("aut")),
     person("Yilong", "Zhang", email = "elong0527@gmail.com", role =  c("aut")),
diff --git a/NAMESPACE b/NAMESPACE
index b2aebbdf..c9b345fe 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -12,9 +12,7 @@ export(mb_weight)
 export(pvalue_maxcombo)
 export(randomize_by_fixed_block)
 export(rpwexp)
-export(rpwexpRcpp)
 export(rpwexp_enroll)
-export(rpwexpinvRcpp)
 export(sim_fixed_n)
 export(sim_pw_surv)
 export(to_sim_pw_surv)
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 8b15b3eb..4e2df0c4 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -1,23 +1,25 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
-#' The piecewise exponential distribution in C++
+#' The piecewise exponential distribution using the inverse CDF method
 #'
 #' @param n Number of observations to be generated.
 #' @param fail_rate A data frame containing `duration` and `rate` variables.
 #'
-#' @export
-rpwexpRcpp <- function(n, fail_rate) {
-    .Call(`_simtrial_rpwexpRcpp`, n, fail_rate)
+#' @noRd
+#'
+rpwexp_inverse_cdf_cpp <- function(n, fail_rate) {
+    .Call(`_simtrial_rpwexp_inverse_cdf_cpp`, n, fail_rate)
 }
 
-#' The piecewise exponential distribution using inverse CDF method in C++
+#' The piecewise exponential distribution using the naive method
 #'
 #' @param n Number of observations to be generated.
 #' @param fail_rate A data frame containing `duration` and `rate` variables.
 #'
-#' @export
-rpwexpinvRcpp <- function(n, fail_rate) {
-    .Call(`_simtrial_rpwexpinvRcpp`, n, fail_rate)
+#' @noRd
+#'
+rpwexp_naive_cpp <- function(n, fail_rate) {
+    .Call(`_simtrial_rpwexp_naive_cpp`, n, fail_rate)
 }
 
diff --git a/R/rpwexp.R b/R/rpwexp.R
index 960b1843..fd5a5932 100644
--- a/R/rpwexp.R
+++ b/R/rpwexp.R
@@ -70,51 +70,5 @@
 rpwexp <- function(
     n = 100,
     fail_rate = data.frame(duration = c(1, 1), rate = c(10, 20))) {
-  n_rate <- nrow(fail_rate)
-
-  if (n_rate == 1) {
-    # Set failure time to Inf if 0 failure rate
-    if (fail_rate$rate == 0) {
-      ans <- rep(Inf, n)
-    } else {
-      # Generate exponential failure time if non-0 failure rate
-      ans <- stats::rexp(n, fail_rate$rate)
-    }
-  } else {
-    # Start of first failure rate interval
-    start_time <- 0
-    # Ends of failure rate interval
-    end_time <- cumsum(fail_rate$duration)
-    # Initiate vector for failure times
-    ans <- rep(0, n)
-    # Index for event times not yet reached
-    indx <- rep(TRUE, n)
-
-    for (i in 1:n_rate) {
-      # Number of event times left to generate
-      n_event_left <- sum(indx)
-
-      # Stop if event is arrived
-      if (n_event_left == 0) {
-        break
-      }
-
-      # Set failure time to Inf for interval i if 0 fail rate
-      if (fail_rate$rate[i] == 0) {
-        ans[indx] <- start_time + rep(Inf, n_event_left)
-      } else {
-        # Generate exponential failure time for interval i if non-0 failure rate
-        ans[indx] <- start_time + stats::rexp(n_event_left, fail_rate$rate[i])
-      }
-
-      # Skip this for last interval as all remaining times are generated there
-      if (i < n_rate) {
-        start_time <- end_time[i]
-        # Update index of event times not yet reached
-        indx <- (ans > end_time[i])
-      }
-    }
-  }
-
-  ans
+  rpwexp_inverse_cdf_cpp(n = n, fail_rate = fail_rate)
 }
diff --git a/R/rpwexp_naive.R b/R/rpwexp_naive.R
new file mode 100644
index 00000000..64b83d3f
--- /dev/null
+++ b/R/rpwexp_naive.R
@@ -0,0 +1,120 @@
+#  Copyright (c) 2023 Merck & Co., Inc., Rahway, NJ, USA and its affiliates.
+#  All rights reserved.
+#
+#  This file is part of the simtrial program.
+#
+#  simtrial is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+#' The piecewise exponential distribution
+#'
+#' The piecewise exponential distribution allows a simple method to specify
+#' a distribution where the hazard rate changes over time.
+#' It is likely to be useful for conditions where failure rates change,
+#' but also for simulations where there may be a delayed treatment effect
+#' or a treatment effect that that is otherwise changing
+#' (for example, decreasing) over time.
+#' `rpwexp_naive()` is to support simulation of both the Lachin and Foulkes (1986)
+#' sample size method for (fixed trial duration) as well as the
+#' Kim and Tsiatis (1990) method (fixed enrollment rates and either
+#' fixed enrollment duration or fixed minimum follow-up);
+#' see [gsDesign::nSurv()].
+#'
+#' Using the `cumulative = TRUE` option, enrollment times that piecewise
+#' constant over time can be generated.
+#'
+#' @param n Number of observations to be generated.
+#' @param fail_rate A data frame containing `duration` and `rate` variables.
+#'   `rate` specifies failure rates during the corresponding interval duration
+#'   specified in `duration`. The final interval is extended to be infinite
+#'   to ensure all observations are generated.
+#'
+#' @noRd
+#'
+#' @examples
+#' # Example 1
+#' # Exponential failure times
+#' x <- rpwexp_naive(
+#'   n = 10000,
+#'   fail_rate = data.frame(rate = 5, duration = 1)
+#' )
+#' plot(sort(x), (10000:1) / 10001,
+#'   log = "y", main = "Exponential simulated survival curve",
+#'   xlab = "Time", ylab = "P{Survival}"
+#' )
+#'
+#' # Example 2
+#'
+#' # Get 10k piecewise exponential failure times.
+#' # Failure rates are 1 for time 0 to 0.5, 3 for time 0.5 to 1, and 10 for > 1.
+#' # Intervals specifies duration of each failure rate interval
+#' # with the final interval running to infinity.
+#' x <- rpwexp_naive(
+#'   n = 1e4,
+#'   fail_rate = data.frame(rate = c(1, 3, 10), duration = c(.5, .5, 1))
+#' )
+#' plot(sort(x), (1e4:1) / 10001,
+#'   log = "y", main = "PW Exponential simulated survival curve",
+#'   xlab = "Time", ylab = "P{Survival}"
+#' )
+rpwexp_naive <- function(
+    n = 100,
+    fail_rate = data.frame(duration = c(1, 1), rate = c(10, 20))) {
+  n_rate <- nrow(fail_rate)
+
+  if (n_rate == 1) {
+    # Set failure time to Inf if 0 failure rate
+    if (fail_rate$rate == 0) {
+      ans <- rep(Inf, n)
+    } else {
+      # Generate exponential failure time if non-0 failure rate
+      ans <- stats::rexp(n, fail_rate$rate)
+    }
+  } else {
+    # Start of first failure rate interval
+    start_time <- 0
+    # Ends of failure rate interval
+    end_time <- cumsum(fail_rate$duration)
+    # Initiate vector for failure times
+    ans <- rep(0, n)
+    # Index for event times not yet reached
+    indx <- rep(TRUE, n)
+
+    for (i in 1:n_rate) {
+      # Number of event times left to generate
+      n_event_left <- sum(indx)
+
+      # Stop if event is arrived
+      if (n_event_left == 0) {
+        break
+      }
+
+      # Set failure time to Inf for interval i if 0 fail rate
+      if (fail_rate$rate[i] == 0) {
+        ans[indx] <- start_time + rep(Inf, n_event_left)
+      } else {
+        # Generate exponential failure time for interval i if non-0 failure rate
+        ans[indx] <- start_time + stats::rexp(n_event_left, fail_rate$rate[i])
+      }
+
+      # Skip this for last interval as all remaining times are generated there
+      if (i < n_rate) {
+        start_time <- end_time[i]
+        # Update index of event times not yet reached
+        indx <- (ans > end_time[i])
+      }
+    }
+  }
+
+  ans
+}
diff --git a/R/sim_pw_surv.R b/R/sim_pw_surv.R
index 70dde971..0c744b2a 100644
--- a/R/sim_pw_surv.R
+++ b/R/sim_pw_surv.R
@@ -174,11 +174,11 @@ sim_pw_surv <- function(
   for (sr in unique_stratum) {
     for (tr in unique_treatment) {
       indx <- x$stratum == sr & x$treatment == tr
-      x$fail_time[indx] <- rpwexpinvRcpp(
+      x$fail_time[indx] <- rpwexp(
         n = sum(indx),
         fail_rate = fail_rate[fail_rate$stratum == sr & fail_rate$treatment == tr, , drop = FALSE]
       )
-      x$dropout_time[indx] <- rpwexpinvRcpp(
+      x$dropout_time[indx] <- rpwexp(
         n = sum(indx),
         fail_rate = dropout_rate[dropout_rate$stratum == sr & dropout_rate$treatment == tr, , drop = FALSE]
       )
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 38308c6e..c7e828f1 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -27,8 +27,6 @@ reference:
 - title: "Distributions"
   contents:
     - rpwexp
-    - rpwexpRcpp
-    - rpwexpinvRcpp
     - rpwexp_enroll
     - fit_pwexp
 
diff --git a/man/rpwexpRcpp.Rd b/man/rpwexpRcpp.Rd
deleted file mode 100644
index 0cee78c8..00000000
--- a/man/rpwexpRcpp.Rd
+++ /dev/null
@@ -1,16 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/RcppExports.R
-\name{rpwexpRcpp}
-\alias{rpwexpRcpp}
-\title{The piecewise exponential distribution in C++}
-\usage{
-rpwexpRcpp(n, fail_rate)
-}
-\arguments{
-\item{n}{Number of observations to be generated.}
-
-\item{fail_rate}{A data frame containing \code{duration} and \code{rate} variables.}
-}
-\description{
-The piecewise exponential distribution in C++
-}
diff --git a/man/rpwexpinvRcpp.Rd b/man/rpwexpinvRcpp.Rd
deleted file mode 100644
index 7335fa0b..00000000
--- a/man/rpwexpinvRcpp.Rd
+++ /dev/null
@@ -1,16 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/RcppExports.R
-\name{rpwexpinvRcpp}
-\alias{rpwexpinvRcpp}
-\title{The piecewise exponential distribution using inverse CDF method in C++}
-\usage{
-rpwexpinvRcpp(n, fail_rate)
-}
-\arguments{
-\item{n}{Number of observations to be generated.}
-
-\item{fail_rate}{A data frame containing \code{duration} and \code{rate} variables.}
-}
-\description{
-The piecewise exponential distribution using inverse CDF method in C++
-}
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index bfcc80ca..64160634 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -10,34 +10,34 @@ Rcpp::Rostream<true>&  Rcpp::Rcout = Rcpp::Rcpp_cout_get();
 Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
 #endif
 
-// rpwexpRcpp
-NumericVector rpwexpRcpp(int n, DataFrame fail_rate);
-RcppExport SEXP _simtrial_rpwexpRcpp(SEXP nSEXP, SEXP fail_rateSEXP) {
+// rpwexp_inverse_cdf_cpp
+NumericVector rpwexp_inverse_cdf_cpp(int n, DataFrame fail_rate);
+RcppExport SEXP _simtrial_rpwexp_inverse_cdf_cpp(SEXP nSEXP, SEXP fail_rateSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< int >::type n(nSEXP);
     Rcpp::traits::input_parameter< DataFrame >::type fail_rate(fail_rateSEXP);
-    rcpp_result_gen = Rcpp::wrap(rpwexpRcpp(n, fail_rate));
+    rcpp_result_gen = Rcpp::wrap(rpwexp_inverse_cdf_cpp(n, fail_rate));
     return rcpp_result_gen;
 END_RCPP
 }
-// rpwexpinvRcpp
-NumericVector rpwexpinvRcpp(int n, DataFrame fail_rate);
-RcppExport SEXP _simtrial_rpwexpinvRcpp(SEXP nSEXP, SEXP fail_rateSEXP) {
+// rpwexp_naive_cpp
+NumericVector rpwexp_naive_cpp(int n, DataFrame fail_rate);
+RcppExport SEXP _simtrial_rpwexp_naive_cpp(SEXP nSEXP, SEXP fail_rateSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< int >::type n(nSEXP);
     Rcpp::traits::input_parameter< DataFrame >::type fail_rate(fail_rateSEXP);
-    rcpp_result_gen = Rcpp::wrap(rpwexpinvRcpp(n, fail_rate));
+    rcpp_result_gen = Rcpp::wrap(rpwexp_naive_cpp(n, fail_rate));
     return rcpp_result_gen;
 END_RCPP
 }
 
 static const R_CallMethodDef CallEntries[] = {
-    {"_simtrial_rpwexpRcpp", (DL_FUNC) &_simtrial_rpwexpRcpp, 2},
-    {"_simtrial_rpwexpinvRcpp", (DL_FUNC) &_simtrial_rpwexpinvRcpp, 2},
+    {"_simtrial_rpwexp_inverse_cdf_cpp", (DL_FUNC) &_simtrial_rpwexp_inverse_cdf_cpp, 2},
+    {"_simtrial_rpwexp_naive_cpp", (DL_FUNC) &_simtrial_rpwexp_naive_cpp, 2},
     {NULL, NULL, 0}
 };
 
diff --git a/src/rpwexpRcpp.cpp b/src/rpwexpRcpp.cpp
deleted file mode 100644
index 6e752aa5..00000000
--- a/src/rpwexpRcpp.cpp
+++ /dev/null
@@ -1,57 +0,0 @@
-#include <Rcpp.h>
-#include <random>
-using namespace Rcpp;
-
-//' The piecewise exponential distribution in C++
-//'
-//' @param n Number of observations to be generated.
-//' @param fail_rate A data frame containing `duration` and `rate` variables.
-//'
-//' @export
-// [[Rcpp::export]]
-NumericVector rpwexpRcpp(int n,
-                      DataFrame fail_rate) {
-  NumericVector duration = fail_rate["duration"];
-  NumericVector rate = fail_rate["rate"];
-  int n_rates = duration.size();
-
-  // Initialize failure times to Inf
-  NumericVector times(n, R_PosInf);
-
-  if (n_rates == 1) {
-    if (rate[0] != 0) {
-      times = rexp(n, rate[0]); // generate exponential failure time if non-0 failure rate
-    }
-  } else {
-    double starttime = 0; // start of first failure rate interval
-    LogicalVector indx(n, 1); // index for event times not yet reached
-    int nindx = n; // number of event times left to generate
-
-    for (int i = 0; i < n_rates; i++) {
-      if (is_false(any(indx))) break; // stop if finished
-      if (rate[i] == 0) {
-        NumericVector temp(nindx, R_PosInf); // set failure time to Inf for inveral i if 0 fail rate
-        times[indx] = temp;
-      } else {
-        NumericVector temp = rexp(nindx, rate[i]); // generate exponential failure time for interval i if non-0 faiurel rate
-        int p = 0;
-        for (int j = 0; j < n; j++) {
-          if (indx[j]) times[j] = temp[p++] + starttime;
-        }
-      }
-
-      if (i < n_rates) { // skip this for last interval as all remaining times are generated there
-        starttime += duration[i]; // update start time for next interval
-        for (int j = 0; j < n; j++) {
-          if (indx[j] == 1 && times[j] <= starttime) {
-            indx[j] = 0; // update index of event times not yet reached
-            nindx--; // update number of event times left to generate
-          }
-        }
-      }
-
-    }
-  }
-
-  return times;
-}
diff --git a/src/rpwexpinvRcpp.cpp b/src/rpwexp_inverse_cdf.cpp
similarity index 69%
rename from src/rpwexpinvRcpp.cpp
rename to src/rpwexp_inverse_cdf.cpp
index 68976719..569a9833 100644
--- a/src/rpwexpinvRcpp.cpp
+++ b/src/rpwexp_inverse_cdf.cpp
@@ -3,16 +3,16 @@
 #include <algorithm>
 using namespace Rcpp;
 
-//' The piecewise exponential distribution using inverse CDF method in C++
+//' The piecewise exponential distribution using the inverse CDF method
 //'
 //' @param n Number of observations to be generated.
 //' @param fail_rate A data frame containing `duration` and `rate` variables.
 //'
-//' @export
+//' @noRd
+//'
 // [[Rcpp::export]]
-NumericVector rpwexpinvRcpp(int n,
-                            DataFrame fail_rate) {
-
+NumericVector rpwexp_inverse_cdf_cpp(int n, DataFrame fail_rate)
+{
   NumericVector duration = fail_rate["duration"];
   NumericVector rate = fail_rate["rate"];
   int n_rates = duration.size();
@@ -23,14 +23,16 @@ NumericVector rpwexpinvRcpp(int n,
   // Get number of piecewise rates
   NumericVector cumTime(n_rates);
   NumericVector cumHaz(n_rates);
-  for (int i = 1; i < n_rates; i++) {
-    cumTime[i] = cumTime[i-1] + duration[i-1];
-    cumHaz[i] = cumHaz[i-1] + duration[i-1] * rate[i-1];
+  for (int i = 1; i < n_rates; i++)
+  {
+    cumTime[i] = cumTime[i - 1] + duration[i - 1];
+    cumHaz[i] = cumHaz[i - 1] + duration[i - 1] * rate[i - 1];
   }
 
   NumericVector::iterator pos;
   int j;
-  for (int i = 0; i < n; i++) {
+  for (int i = 0; i < n; i++)
+  {
     pos = std::upper_bound(cumHaz.begin(), cumHaz.end(), times[i]);
     j = pos - cumHaz.begin() - 1;
     times[i] = cumTime[j] + (times[i] - cumHaz[j]) / rate[j];
diff --git a/src/rpwexp_naive.cpp b/src/rpwexp_naive.cpp
new file mode 100644
index 00000000..de73a841
--- /dev/null
+++ b/src/rpwexp_naive.cpp
@@ -0,0 +1,71 @@
+#include <Rcpp.h>
+#include <random>
+using namespace Rcpp;
+
+//' The piecewise exponential distribution using the naive method
+//'
+//' @param n Number of observations to be generated.
+//' @param fail_rate A data frame containing `duration` and `rate` variables.
+//'
+//' @noRd
+//'
+// [[Rcpp::export]]
+NumericVector rpwexp_naive_cpp(int n, DataFrame fail_rate)
+{
+  NumericVector duration = fail_rate["duration"];
+  NumericVector rate = fail_rate["rate"];
+  int n_rates = duration.size();
+
+  // Initialize failure times to Inf
+  NumericVector times(n, R_PosInf);
+
+  if (n_rates == 1)
+  {
+    if (rate[0] != 0)
+    {
+      times = rexp(n, rate[0]); // Generate exponential failure time if non-zero failure rate
+    }
+  }
+  else
+  {
+    double starttime = 0;     // Start of first failure rate interval
+    LogicalVector indx(n, 1); // Index for event times not yet reached
+    int nindx = n;            // Number of event times left to generate
+
+    for (int i = 0; i < n_rates; i++)
+    {
+      if (is_false(any(indx)))
+        break; // stop if finished
+      if (rate[i] == 0)
+      {
+        NumericVector temp(nindx, R_PosInf); // Set failure time to Inf for interval i if 0 fail rate
+        times[indx] = temp;
+      }
+      else
+      {
+        NumericVector temp = rexp(nindx, rate[i]); // Generate exponential failure time for interval i if non-0 failure rate
+        int p = 0;
+        for (int j = 0; j < n; j++)
+        {
+          if (indx[j])
+            times[j] = temp[p++] + starttime;
+        }
+      }
+
+      if (i < n_rates)
+      {                           // Skip this for last interval as all remaining times are generated there
+        starttime += duration[i]; // Update start time for next interval
+        for (int j = 0; j < n; j++)
+        {
+          if (indx[j] == 1 && times[j] <= starttime)
+          {
+            indx[j] = 0; // Update index of event times not yet reached
+            nindx--;     // Update number of event times left to generate
+          }
+        }
+      }
+    }
+  }
+
+  return times;
+}
diff --git a/tests/testthat/test-independent_test_rpwexpinvRcpp.R b/tests/testthat/test-independent_test_rpwexp_inverse_cdf_cpp.R
similarity index 62%
rename from tests/testthat/test-independent_test_rpwexpinvRcpp.R
rename to tests/testthat/test-independent_test_rpwexp_inverse_cdf_cpp.R
index 2f1879dd..0bd56680 100644
--- a/tests/testthat/test-independent_test_rpwexpinvRcpp.R
+++ b/tests/testthat/test-independent_test_rpwexp_inverse_cdf_cpp.R
@@ -1,6 +1,6 @@
-test_that("rpwexpinvRcpp handles 0 fail_rate for final period", {
+test_that("rpwexp_inverse_cdf_cpp handles 0 fail_rate for final period", {
   # 0 failure rate for last period
-  s <- simtrial::rpwexpinvRcpp(
+  s <- simtrial:::rpwexp_inverse_cdf_cpp(
     n = 100,
     fail_rate = data.frame(duration = c(0, 2), rate = c(5, 0))
   )
@@ -8,9 +8,9 @@ test_that("rpwexpinvRcpp handles 0 fail_rate for final period", {
   testthat::expect_equal(min(s), Inf)
 })
 
-testthat::test_that("rpwexpinvRcpp handles 0 fail rate properly for one period", {
+testthat::test_that("rpwexp_inverse_cdf_cpp handles 0 fail rate properly for one period", {
   # 0 failure rate
-  s <- simtrial::rpwexpinvRcpp(
+  s <- simtrial:::rpwexp_inverse_cdf_cpp(
     n = 100,
     fail_rate = data.frame(duration = c(1), rate = c(0))
   )
@@ -19,9 +19,9 @@ testthat::test_that("rpwexpinvRcpp handles 0 fail rate properly for one period",
 })
 
 
-testthat::test_that("rpwexpinvRcpp handles 0 fail rate properly for multiple periods", {
+testthat::test_that("rpwexp_inverse_cdf_cpp handles 0 fail rate properly for multiple periods", {
   # 0 failure rate for 1st period (with duration of 1 time unit)
-  s <- simtrial::rpwexpinvRcpp(
+  s <- simtrial:::rpwexp_inverse_cdf_cpp(
     n = 100,
     fail_rate = data.frame(duration = c(1, 2), rate = c(0, 5))
   )