Without posting equations and examples, those who are familiar with statistics should know it always involves differentiating and integrating pdfs, from calculating moments to estimators. Pdf multiplication (joint pdfs, conditional pdfs) results in some complex formulas of several variables. Be it moment-generating-functions, MLE, or Bayes estimators, we integrate the unknown parameter over its domain while other parameters remain undefined. For a simple example, the Bayes estimator for the p parameter of binomial distributions will lead to integrating the conditional pdf involving p, n (sample size) and x(number of success within sample), integrating p while n and x remain unknown constants. This kind of integration is symbolic computation, right? Is there any way to approach it using numerical integration? Point being, playing with R/Matlab/Octave (involving sympy), it seems symbolic computation is a difficult and computationally-inefficient thing to do. A relatively simple textbook example involving integrating a conditional pdf which can be done by hand, took Matlab over 10-mins to complete and I couldn't even tell if it would succeed. I changed something to see how it would affect the result, it took another 15-mins. It isn't feasible. How are people doing these computations in practice?