gamma_sum.R

generate_gamma_sum <- function(N, k, theta) {
	n_means = length(k)
  sums = numeric(N)
  for(n in 1:N) sums[n] = sum(rgamma(n_means, shape = k, scale = theta))
  list(sums = sums, shape = k, scale = theta)
}

k=c(3,4,5)
theta=c(1, 2, 1)
	generate_gamma_sum(10000, c(3,4,5), c(1, 2, 1)) %$% 
	sums %>%
	MASS::fitdistr("gamma") 
	
	k_sum = sum(theta*k)^2/sum(theta^2*k)
	theta_sum = sum(k*theta)/k_sum