Code does simulation 10^6 times

[Alhejaili]

This code does simulation 10^6 times to calculate the CDF for Poisson bivariate gamma distribution by generating random variables from previous distributions. Now I want to apply this code to real data in an excel file. How I can do that?
Explanation of the code: at the beginning of the loop will generate one value from Poisson if this value is zero so go to the next, if not zero then in (xx) generate the values from gamma equal to the value of Poisson ( if Poisson equal 2 then generate 2 value from gamma if Poisson is 4 then generate 4 value from gamma and so on), also in (yy) generate the values from gamma equal to the value of Poisson as explained earlier. But in real data, I want to take the value from the vector that follows Poisson and take the value from other vectors following gamma based on the taken Poisson value as explained above.

s <- c(2.3, 1.9)
l <- 4
a <- c(1, 1)
a1=a[1]
a2=a[2]
b <- c(2, 4)
b1=b[1]
b2=b[2]
count=0
ll=0

for (j in 1:10000000) 
{
  f = rpois(1,l)
  
  if(f == 0)
  {
    count = count + 1
    next
  }
  
  xx = rgamma(f,a1,rate=1/b1)
  yy = rgamma(f,a2,rate=1/b2)
  
  if(sum(xx)== s[1] && sum(yy)<= s[2])
  {
    ll = ll + 1
  }

  if(sum(xx)< s[1] && sum(yy)== s[2])
  {
    ll = ll + 1
  }
  
  if(sum(xx)< s[1] && sum(yy)< s[2])
  {
    count=count+1
  }
  
}

pp1=((count)/10000000)

pp2=((count+ll)/10000000)

pp =((count+0.5*ll)/10000000)

[AlexisW]

It is extremely unclear to me what you're trying to achieve, and whether it's correct, noting that it's essentially impossible for ll to be more than 0 (it might happen exceptionally due to rounding error that ll = 1), but ou can put your code in a function and call it with whatever value of f you want:

s <- c(2.3, 1.9)
l <- 4
a <- c(1, 1)
a1=a[1]
a2=a[2]
b <- c(2, 4)
b1=b[1]
b2=b[2]


do_my_thing <- function(f){
  if(f == 0)
  {
    return(data.frame(count = 1, ll = 0))
  }
  count <- 0
  ll <- 0
  xx = rgamma(f,a1,rate=1/b1)
  yy = rgamma(f,a2,rate=1/b2)
  
  if(sum(xx)== s[1] && sum(yy)<= s[2])
  {
    ll = ll + 1
  }
  
  if(sum(xx)< s[1] && sum(yy)== s[2])
  {
    ll = ll + 1
  }
  
  if(sum(xx)< s[1] && sum(yy)< s[2])
  {
    count <- 1
  }
  return(data.frame(count = count, ll = ll))
}

result <- lapply(rpois(1e4, l), do_my_thing) |>
  do.call(rbind, args = _)

final_ll <- sum(result$ll)
final_count <- sum(result$count)

pp1=((final_count)/10000000)

pp2=((final_count+final_ll)/10000000)

pp =((final_count+0.5*final_ll)/10000000)

Some explanation: I put the content of the loop in a function, at each call of the function we start count and ll at 0 and increment them when needed, and return them as a dataframe. That way the generation of f happens outside the function.

I would note that, since xx and yy are results of rgamma, their sum is extremely unlikely to be exactly equal to an exact value, specifically the probability that sum(xx) == 2.3 is 0 (but since we're working with a computer, there are always rounding errors, so there is a small but non-zero chance that sum(xx) ~= 2.3 within the margin of error). Same for yy. So I'm not sure the computation of ll serves any purpose.

So I think you could further simplify your code, for example this should give you the same result:

get_count <- function(f){
  f == 0 || (sum(rgamma(f,a1,rate=1/b1))< s[1] &&
               sum(rgamma(f,a2,rate=1/b2))< s[2])
}

res <- sapply(rpois(1e4, l), get_count)
sum(res)

where sum(res) is equivalent to count in your for loop.