Normality of a ratio of random variables For the context, I am measuring the execution time of a piece of code. I run it $N$ times with different problem sizes, its duration is proportional to the problem size. I am interested in the performance, in terms of operations per second, and in particular the "average" performance.
Here are my variables:

*

*$F_i$ is the number of operations performed by the run $i$. No assumption is made on these.

*$T_i \sim \alpha F_i + \beta + \epsilon$ is the duration of run $i$, which depends on the number of operations. No assumption is made on $\epsilon$ (it could be normally distributed, or uniformly, or whatever).

*$P_i = F_i/T_i$ is the performance of run $i$, in operations per second.

I have two alternative definitions for the "average" performance. These definitions are not equivalent, the second one makes more sense in my context.

*

*$M_1 = \frac{1}{N}\sum_i P_i$, the arithmetic mean of the individual performance values.

*$M_2 = \sum F_i / \sum T_i$, the ratio of the total number of operations by the total duration.

Now I am interested in the distribution of $M_1$ and $M_2$ if I were to repeat this experiment several times. I think we can expect $M_1$ to be normally distributed, with the central limit theorem. My question is about $M_2$, do we know something about its distribution?

I tried to simulate this scenario with some dummy data. I got this plot. The red plots are the theoretical densities of the normal distribution with the sample mean and standard deviation, the black plots are the real density plots. To me it seems pretty clear that both variables ($M_1$ and $M_2$) follow a normal distribution, but I cannot convince myself why it is the case for $M_2$.

The code:
library(dplyr)
library(ggplot2)
library(patchwork)

N = 1000           # number of calls in an experiment
nb_exp = 10000     # number of experiments
max_flop = 1e10    # maximal number of operations in a call
mean_rate = 7e-11  # inverse of the flop rate
intercept = 1e-2   # minimum duration of a call
small_calls = rep(c(1, 4, 16), 50) # each experiment also has these small calls
set.seed(42)

df = data.frame()
for(i in (1:nb_exp)) {
    exp = data.frame(F=c(runif(N, 1, max_flop), small_calls)) %>%
        mutate(T=rnorm(n(), F*mean_rate, F*mean_rate/20) + intercept) %>%
        mutate(P=F/T)
    mean1 = exp %>% pull(P) %>% mean()
    mean2 = exp %>% summarise(F=sum(F), T=sum(T)) %>% mutate(P=F/T) %>% pull(P)
    df = rbind(df, data.frame(mean1=mean1, mean2=mean2))
}

plot_density <- function(df, col) {
    mu = mean(df[[col]])
    sigma = sd(df[[col]])
    ggplot(df) +
        aes_string(x=col) +
        geom_density() +
        stat_function(fun=dnorm, n=1001, args=list(mean=mu, sd=sigma), color="red", linetype="dashed") +
        theme_void()
}

plot_bivariate <- function(df, x_col, y_col) {
    plot_top = plot_density(df, x_col)

    plot_right = plot_density(df, y_col) +
        coord_flip()

    scatter_plot = ggplot(df) +
        aes_string(x=x_col, y=y_col) +
        geom_point() +
        stat_ellipse() +
        theme_minimal()

    return(
        plot_top + plot_spacer() + scatter_plot + plot_right +
        plot_layout(widths = c(4, 1), heights = c(1, 4)) &
        theme(legend.position='none')
    )
}

plot_bivariate(df, "mean1", "mean2") &
    xlab("Arithmetic mean of individual performances") &
    ylab("Ratio of the total number of operations by the total duration")

 A: Since $\epsilon$ is random, I am assuming it can be different each run and I will put a subscript $i$ on it.
Let $U=\frac{\sum_{i=1}^N{F_i}}{N}$ and $V=\frac{\sum_{i=1}^N{T_i}}{N}$.
The asymptotic variance of $U$ is $\frac{Var(F_1)}{N}$ and the asymptotic variance of $V$ is $\frac{Var(T_1)}{N}=\frac{\alpha^2Var(F_1)+Var(\epsilon_1)}{N}$.
Both $U$ and $V$ are asymptotically normally distributed by the Central Limit Theorem.
Since all variables with different subscripts are independent and have covariance 0, the covariance between the two is
$$Cov(U,V)=Cov \left(\frac{\sum_{i=1}^N{F_i}}{N},\frac{\sum_{i=1}^N{T_i}}{N} \right)=\frac{\sum_{i=1}^N{Cov(T_i,F_i)}}{N^2}$$
$$=\frac{\sum_{i=1}^N{Cov(\alpha F_i+\beta+\epsilon_i,F_i)}}{N^2}=\frac{\alpha \sum_{i=1}^N{ Cov( F_i,F_i)}}{N^2}=\frac{\alpha \sum_{i=1}^N{ Var(F_i)}}{N^2}=\frac{\alpha Var(F_1)}{N}$$
To summarize the above
$$\sqrt{N}\left(\begin{bmatrix}U\\V\end{bmatrix}-
\begin{bmatrix}E [F_1]\\\alpha E [F_1]+\beta+E[\epsilon_1]\end{bmatrix}
\right)
\rightarrow N\left(\begin{bmatrix}0\\0\end{bmatrix},\begin{bmatrix}Var(F_1) & \alpha Var(F_1)\\\alpha Var(F_1) & \alpha^2Var(F_1)+Var(\epsilon_1)\end{bmatrix}\right)$$
Since $M_2=\frac{U}V$, using the multivariate Delta method, we have
$$\sqrt{N}\left(M_2-\frac{E [F_1]}{\alpha E [F_1]+\beta+E[\epsilon_1]}\right)
\rightarrow 
N\left(0,
g'
\begin{bmatrix}Var(F_1) & \alpha Var(F_1)\\\alpha Var(F_1) & \alpha^2Var(F_1)+Var(\epsilon_1)\end{bmatrix} g
\right)$$
where $g$ is the gradient vector $\begin{bmatrix}\frac{1}{\alpha E [F_1]+\beta+E[\epsilon_1]}\\-\frac{E [F_1]}{(\alpha E [F_1]+\beta+E[\epsilon_1])^2}\end{bmatrix}$
