This paper provides a numerical framework for computing the achievable rate region of memoryless multiple access channel (MAC) with a continuous alphabet from data. In particular, we use recent results on variational lower bounds on mutual information and KL-divergence to compute the boundaries of the rate region of MAC using a set of functions parameterized by neural networks. Our method relies on a variational lower bound on KL-divergence and an upper bound on KL-divergence based on the f-divergence inequalities. Unlike previous work, which computes an estimate on mutual information, which is neither a lower nor an upper bound, our method estimates a lower bound on mutual information. Our numerical results show that the proposed method provides tighter estimates compared to the MINE-based estimator at large SNRs while being computationally more efficient. Finally, we apply the proposed method to the optical intensity MAC and obtain a new achievable rate boundary tighter than prior works.