Summations zu.mpfr and zstar based on a probability sample

Question

Here is my code:

library(gtools)
library(Rmpfr)
library(OBsMD)
precBits=1000000
K=9
n=1000
p=mpfr(c(K:1),precBits = precBits);p=mpfr(p,precBits = precBits)/mpfr(sum(p),precBits = precBits)
#p is a vector of probabilities

zu.mpfr <- function(freq, u, precBits){
  sum=mpfr(0,precBits = precBits)
  n=sum(freq)
  for (i in (1:length(freq))){
    prod=mpfr(1,precBits = precBits)
    for (j in (0:(u-1))){
      prod=prod*(mpfr(freq[i]-j,precBits = precBits))/mpfr(n-j,precBits=precBits)
    }
    sum=sum+prod
  }
  return(sum)
}

etastar<-rep(mpfr(0,precBits = precBits),K-1)
etastar=c(etastar)
for (i in 1:(K-1)){
  etastar[i]=sum(p^mpfr(i+1,precBits))
}

In the above codes, p is a vector of probabilities. zu.mpfr is \$z_u =\sum_k \prod_{j=0}^{u-1}\frac{Y_k-j}{n-j},\$ where \$u\$ can be any positive integer (but no more than \$n\$), and \$Y_k\$ is the observed frequency in a random sample for catagory \$k\$. The inputs for zu.mpfr are the sample frequency counts, \$u\$, and precBits (the preset accuracy). For the frequency, for example, if a sample is \$\{1,1,2,3,3,3,4\}\$, then the corresponding frequencies are \$\{2,1,3,1\}\$, since observation "1" appears twice, "2" appears once, "3" appears three times, and "4" appears once. The frequency is directly generated by rmultinom in the next paragraph of codes. etastar is a vector of length \$K-1\$. And the \$i\$-th element, \$\eta^{*}(i)\$, is \$\sum_k p_k ^{i+1}.\$ As a summary, zu.mpfr is a value calculated from a random sample, and etastar is a vector calculated from the probability distribution p.

Following is the part that I need to make it faster:

start.time<-Sys.time()
sample=rmultinom(1,n,asNumeric(p))
zstar<-rep(mpfr(0,precBits = precBits),K-1)
zstar=c(zstar)
for(i in 1:(K-1)){
  zstar[i]=zu.mpfr(sample,i+1,precBits)
}

end_time <- Sys.time()
end_time-start.time

In the latter codes, I used rmultinom to get the sample, which is the frequency counts from a random sample that generated by the probability distribution p. Then I calculated zstar, which is a vector of length \$K-1\$, with the \$i\$-th element zu.mpfr(sample, i+1, precBits). Here comes my problem. It takes too long to calculate the vector zstar. By timing it, "Time difference of 6.207007 secs" on my end. I am wondering if I can make it better. I need it faster because I need to run the simulation with thousands of iterations. I am not sure if the slowness is caused by mpfr but I do need it to enforce accuracy.

Please let me know if any further clarification is needed.

Answer 1

First suggestion would be to write a function for the high precision conversion because you keep using mpfr(..., precBits = precBits) over and over:

Z <- function(x) mpfr(x, precBits = 1000000)

It won't make your code faster but it will save you a lot of typing and make it easier on the eye.

Next, if you do not know why your code is slow, you need to learn to use the profiler:

Rprof(tmp <- tempfile())
[INSERT_YOUR_CODE_HERE]
Rprof()
summaryRprof(tmp)
unlink(tmp)

It will show you that 90% of the time is spent doing the division in your prod = prod * ... line of code. So a solution would be to compute the numerator (num) and denominator (den) by iteration, then do a single division at the end:

zu.mpfr <- function(freq, u) {
  tot <- Z(0)
  n <- sum(freq)
  for (i in seq_along(freq)) {
    num <- Z(1)
    den <- Z(1)
    for (j in (0:(u-1))) {
      num <- num * Z(freq[i] - j)
      den <- den * Z(n - j)
    }
    prod <- num / den
    tot <- tot + prod
  }
  return(tot)
}

This alone makes your code three to four times faster.

Next, I would recommend looking for places where you can vectorize. I understand that your code is dealing with high precision numbers so maybe there is no other way than using for loops if only the simple operators (+, *, /) are available... However, it appears that if x is a vector of mpfr numbers, prod(x) also returns a mpfr number so I am ready to trust it is doing the right thing. If so, then you could do:

zu.mpfr <- function(freq, u) {
  tot <- Z(0)
  n <- sum(freq)
  J <- 0:(u-1)
  for (i in seq_along(freq)) {
    tot <- tot + prod(Z(freq[i] - J)) / prod(Z(n - J))
  }
  return(tot)
}

Have you noticed that we are computing the same denominator over and over? Move it outside the loop to save a little more:

zu.mpfr <- function(freq, u) {
  tot <- Z(0)
  n <- sum(freq)
  J <- 0:(u-1)
  den <- prod(Z(n - J))
  for (i in seq_along(freq)) {
    tot <- tot + prod(Z(freq[i] - J)) / den
  }
  return(tot)
}

Last, if memory size is not an issue (i.e. K remains small), you could save even more time by working with matrices. It's a bit harder to digest (including the computation of row-wise cumulative products), but it's elegant how short your code becomes (notice we don't need zu.mpfr anymore). This version is 5 to 6 times faster than your original one:

  sample <- c(rmultinom(1, n, asNumeric(p)))
  x <- Z(outer(sample, 0:(K-1), FUN = "-"))
  y <- Z(sum(sample) - 0:(K-1))
  z <- apply(x, 1, cumprod) / cumprod(y)
  zstar <- tail(rowSums(z), -1)

Answer 2

After some testing, I would suggest that you decrease the precBits, because that changes a lot:

precBits = 10000

I did not see any (significant) changes in results using 10000.

Also, you can try this function, which uses one less loop and could be faster for larger K values:

FUN2 <- function(freq, u, precBits) {
  sum = mpfr(0, precBits = precBits)
  n = sum(freq)
  prod = mpfr(1, precBits = precBits)
  j <- (0:(u - 1))
  for (i in (1:length(freq))) {
      v1 <- mpfr(freq[i] - j, precBits = precBits)
      v2 <- mpfr(n - j, precBits = precBits)
      prod = prod(v1 / v2)
      sum = sum + prod
    }
  return(sum)
}
zstar2 <- rep(mpfr(0, precBits = precBits), K - 1)
for (i in 1:(K - 1)) {
  zstar2[i] = FUN2(sample, i + 1, precBits)
}
all.equal(zstar, zstar2) == T
zstar - zstar2

or this without loops:

FUN4 <- function(freq, u, precBits) {
  n = sum(freq)
  j <- (0:(u - 1))
  xx <- t(apply(freq, 1, function(x) x - j))
  v1 <- mpfr(xx, precBits = precBits)
  v2 <- mpfr(n - j, precBits = precBits)
  xx <- v1 / v2
  prods <- apply(xx, 1, prod)
  s <- sum(prods)
  return(s)
}

I suggest that you try them out with different precBits values and compare the results.