Simulate vectors conditional on custom distribution

r, simulation

Solution

As far as I understand what you're after, I would start by converting `den` to an `rle` object. (Here using data from your Example 3)

EDIT: Add 100% at day 364 to `den`

if(max(den)!=364) den <- c(den, 364)
(rleDen <- rle(den))
# Run Length Encoding
#   lengths: int [1:92] 3 1 1 1 1 1 1 1 1 1 ...    # 92 intervals
#   values : num [1:92] 1 2 4 6 8 9 12 15 17 21 ...
percDur <- rleDen$lengths            # Percentage of total duration in each interval
atDay <- rleDen$values               # What day that percentage was reached
intWidth <- diff(c(0, atDay), k = 1) # Interval width
durPerDay <- 1440                    # Max observation time per day
percPerDay <- durPerDay/a*100        # Max percentage per day
cumPercDur <- cumsum(percDur)        # Cumulative percentage in each interval
maxPerInt <- pmin(percPerDay * diff(c(0, atDay), 1),
  percDur + 1)                       # Max percent observation per interval

set.seed(1)
nsims <- 10                          # Desired number of simulations
sampMat <- matrix(0, ncol = length(percDur), nrow = nsims) # Matrix to hold sim results

To allow for randomness while considering the limitation of a maximum 1440 minutes of observation per day, check to see if there are any long intervals (i.e., any intervals in which the jump in percentage cannot be completely achieved in that interval)

if(any(percDur > maxPerInt)){
  longDays <- percDur > maxPerInt
  morePerInt <- maxPerInt - percDur
  perEnd <- c(which(diff(longDays,1) < 0), length(longDays))
# Group intervals into periods bounded by "long" days
# and determine if there are any long periods (i.e., where
# the jump in percentage can't be achieved in that period)
  perInd <- rep(seq_along(perEnd), diff(c(0, perEnd)))
  perSums <- tapply(percDur, perInd, sum)
  maxPerPer <- tapply(maxPerInt, perInd, sum)
  longPers <- perSums > maxPerPer
# If there are long periods, determine, starting with the last period, when the
# excess can be covered. Each group of periods is recorded in the persToWatch
# object
  if(any(longPers)) {
    maxLongPer <- perEnd[max(which(longPers))]
    persToWatch <- rep(NA, length(maxLongPer))
    for(kk in rev(seq_len(maxLongPer))) {
      if(kk < maxLongPer && min(persToWatch, na.rm = TRUE) <= kk) next
        theSums <- cumsum(morePerInt[order(seq_len(kk),
          decreasing = TRUE)])
        above0 <- which(rev(theSums) > 0)
        persToWatch[kk] <- max(above0[which(!perInd[above0] %in% c(perInd[kk],
          which(longPers)) & !above0 %in% which(longDays))])
    }
  }
}

Now we can start the randomness. The first component of the sampling determines the overall proportion of `a` that occurs in each of the intervals. How much? Let `runif` decide. The upper and lower limits must reflect the maximum observation time per day and the excess amount of any long days and periods

  for(jj in seq_along(percDur[-1])) {
    upperBound <- pmin(sampMat[, jj] + maxPerInt[jj],
      cumPercDur[jj] + 1)
    lowerBound <- cumPercDur[jj]
# If there are long days, determine the interval over which the
# excess observation time may be spread
    if(any(percDur > maxPerInt) && any(which(longDays) >= jj)) {
      curLongDay <- max(which(perInd %in% perInd[jj]))
      prevLongDay <- max(0, min(which(!longDays)[which(!longDays) <= jj]))
      curInt <- prevLongDay : curLongDay
# If there are also long periods, determine how much excess observation time there is
      if(any(longPers) && maxLongPer >= jj) {
        curLongPerHigh <- min(which(!is.na(persToWatch))[
          which(!is.na(persToWatch)) >= jj])
        curLongPerLow <- persToWatch[curLongPerHigh]
        longInt <- curLongPerLow : curLongPerHigh
        curExtra <- max(0,
          cumPercDur[curLongPerHigh] - 
          sum(maxPerInt[longInt[longInt > jj]]) - 
          sampMat[, jj, drop = FALSE])
      } else {
        curExtra <- cumPercDur[curLongDay] - 
          (sum(maxPerInt[curInt[curInt > jj]]) +
          sampMat[, jj, drop = FALSE])
      }
# Set the lower limit for runif appropriately
      lowerBound <- sampMat[, jj, drop = FALSE] + curExtra
    }
# There may be tolerance errors when the observations are tightly
# packed
    if(any(lowerBound - upperBound > 0)) { 
      if(all((lowerBound - upperBound) <= .Machine$double.eps*2*32)) {
        upperBound <- pmax(lowerBound, upperBound)
      } else {
        stop("\nUpper and lower bounds are on the wrong side of each other\n",
          jj,max(lowerBound - upperBound))
      }
    }
    sampMat[, jj + 1] <- runif(nsims, lowerBound, upperBound)
  }

Then add 100 percent to the end of the results and calculate the interval-specific percentage

  sampMat2 <- cbind(sampMat[, seq_along(percDur)], 100)
  sampPercDiff <- t(apply(sampMat2, 1, diff, k = 1))

The second component of the randomness determines the distribution of `sampPercDiff` over the interval widths `intWidth`. This still requires more thought in my opinion. For instance, how long does a typical episode last compared to the unit of time under consideration?

For each interval, determine if the random percentage needs to be allocated over multiple time units (in this case days). EDIT: Changed the following code to limit percentage increase when `intWidth > 1`.

library(foreach)
  ep.dur<-foreach(ii = seq_along(intWidth),.combine=cbind)%do%{
    if(intWidth[ii]==1){
      ret<-sampPercDiff[, ii, drop = FALSE] * a / 100
      dimnames(ret)<-list(NULL,atDay[ii])
      ret
    } else {
      theDist<-matrix(numeric(0), ncol = intWidth[ii], nrow = nsims)
      for(jj in seq_len(intWidth[ii]-1)){
        theDist[, jj] <- floor(runif(nsims, 0, pmax(0,
          min(sampPercDiff[, ii], floor(sampMat2[,ii + 1])-.Machine$double.eps -
          sampMat2[,ii]) * a / 100 - rowSums(theDist, na.rm = TRUE))))
      }
      theDist[, intWidth[ii]] <- sampPercDiff[, ii] * a / 100 - rowSums(theDist,
        na.rm = TRUE)
      distOrder <- replicate(nsims, c(sample.int(intWidth[ii] - 1),
        intWidth[ii]), simplify = FALSE)
      ret <- lapply(seq_len(nrow(theDist)), function(x) {
        theDist[x, order(distOrder[[x]])]
      })
      ans <- do.call(rbind, ret)
      dimnames(ans) <- list(NULL, atDay[ii]-((intWidth[ii]:1)-1))
      ans
    }
  }

The duration time is sampled randomly for each time unit (day) in the interval to which it is to be distributed. After breaking up the total duration into daily observed times, these are then assigned randomly to the days in the interval.

Then, multiply the sampled and distributed percentages by `a` and divide by 100

ep.dur[1, 1 : 6]
#         1         2         3         4         5         6 
# 1095.4475  315.4887    1.0000  578.9200   13.0000  170.6224 

ncol(ep.dur)
# [1] 364

apply(ep.dur, 1, function(x) length(which(x == 0)))
# [1] 131 133 132 117 127 116 139 124 124 129

rowSums(ep.dur)/a
# [1] 1 1 1 1 1 1 1 1 1 1

plot(ep.dur[1, ], type = "h", ylab = "obs time")

Problem

Original source