# Code does simulation 10^6 times

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)

``````

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.

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.