1

I've been trying to create a function for the Irwin Hall distribution that doesn't face the same issue as the unifed package implementation. Because the function suffers from numerical issues, I figured using Rmpfr to be able to specify precision would help. I was wrong. I also attempted to use as.brob() from the Brobdingnag package, but it was very slow (I gave up waiting for it).

Below is my attempt at the solution using the formula from the unifed package as explained in this glorious answer to the sum of uniform random variables question.

library("unifed")
library("Rmpfr")

bits<-260
dirwin.hall2 <- function(x, n, a=0, b=1){


x is a mean from unif(a,b)


we want it to look like a sum from unif(0,1)


scaled.x <- xn               # sum instead of mean
  scaled.x <- scaled.x - na    # shift back to unif(0,b-a)
  scaled.x <- scaled.x/(b-a)    # scale back to unif(0,1)


f <- rep(NA, length(x))
  for(i in 1:length(x)){
    ret<- mpfr(0,bits)
    for (k in 0:floor(scaled.x[i])){ 
     ret <- ret + mpfr((-1)^k, bits) * mpfr(choose(n, k),bits) * mpfr((scaled.x[i] - k)^(n - 1),bits)
    }
    f[i]<-ret
  }
  suppressWarnings((n.samp/(b-a))(1/factorial(n-1))  sapply(f, asNumeric))
}


dirwin.hall(35,50)
dirwin.hall(36,50)
dirwin.hall(37,50)
dirwin.hall(38,50)


dirwin.hall2(35/50,50)
dirwin.hall2(36/50,50)
dirwin.hall2(37/50,50)
dirwin.hall2(38/50,50)
user1329307
  • 113
  • Which “the same issue”? What exactly are you trying to do? You only told us that you are to implement some algorithm to solve some issue, but without saying what are they. – Tim Mar 08 '23 at 21:49

0 Answers0