Numerically finding confidence interval bounds I am asking an R question here on the basis that statistical expertise is needed, i.e., 

If the language is statistically oriented (such as R, SAS, Stata,
  SPSS, etc.), then decide based on the nature of your question: if it
  needs statistical expertise to understand or answer, ask it here

I have a symmetric (about some value $\mu$) probability density function $f$ with support in $(-\infty, \infty)$. Let's suppose, for simplicity's sake, that it's the normal distribution $$f(x; \mu, \sigma) = \dfrac{1}{\sigma\sqrt{2\pi}}\exp\left[-\left(\frac{x-\mu}{\sigma}\right)^2/2\right]\text{.}$$
Ignoring the fact that there's a dnorm function in R, let's say I wrote this out as a function:
f <- function(x, mu, sigma){
1/(sigma*sqrt(2*pi))*exp(-((x-mu)/sigma)^2/2)
}

I would like to solve the equation
$$\int_{a}^{b}f(x; \mu, \sigma)\text{ d}x = \pi \in [0, 1]\text{ fixed.}$$
where $\mu = \dfrac{a+b}{2}$, for $a$ and $b$. 
Without standardizing the normal distribution and without using any of the "normal"-supplied functions in R (e.g., dnorm, pnorm, qnorm, rnorm and related functions), how do I solve for $a$ and $b$ in R? 
Here's what I do know:


*

*R has an integrating function. But I'm not sure how to use this to help me solve for the bounds.

*The normal distribution PDF, along with the PDF of distribution I am working with, do not have closed forms for the CDF.


I would guess that simulation is involved, but I wouldn't know where to start.
(Note: I know this seems silly, but I would like to know how this works for the normal distribution so that I can extend it to a general PDF.)
 A: Yes use integrate to evaluate the integral and then e.g. uniroot to solve $\int_{\mu-a}^{\mu+a} f(x)dx - \pi = 0$ as follows (n.b. I changed the parameterisation slightly):
f <- function(x, mu, sigma) {
    1/(sigma*sqrt(2*pi)) * exp(-((x-mu)/sigma)^2/2)
}

g<-function(a,mu,sigma,p) {
    integrate(function(x){f(x,mu,sigma)}, mu-a, mu+a)$value - p
}

uniroot(g, interval=c(-10,10), mu=0,sigma=1, p =0.95)$root

Output:
[1] 1.959969

Note that the (-10,10) is the search interval.  Check ?uniroot for more info.
N.b. no need for simulation.  Also, even if a closed form solution is not available, many distributions have good numerical approximations to their inverses.
A: I am not familiar with R. 
But a (simple) solution is based on Monte Carlo:
You can draw n samples from the distribution. 
With the following algorithm you can find a and b:  


*

*Divide the samples in two equally sized arrays: 

*the samples greater as the median (gm) and 

*the samples lower than the median (lm).   

*to get a : sort (lm) descent. a is approximate the value of the $n * \pi/2$-item ($\pi$ is the value of the total integral)

*analog for b: sorting order asc

