How to simulate random observations from a specific distribution? I am asking for a general approach about how to construct algorithms akin to, for example, the rnorm function in R given that one has, say, a closed-form probability distribution function and corresponding cumulative distribution function (and quantile function, too, if needed).
How to go about it?
I am not interested in the normal distribution per se, nor interested in a specific language, but the approach.
 A: The Art Of Computer Programming Vol.2 Chapter 3 has all the different methods used. I highly recommend reading the entire chapter.
The simplest and most obvious method is by using the inverse cumulative distribution function (CDF): $F^{-1}(p)$, which is also the quantile function. Generate uniform random numbers and plug them into the inverse CDF, you'll get what you're looking for. This is what I usually do in Excel NORMS.INV(RAND()), because I'm too lazy to code Box-Mueller myself.
This is not the method used for Normal distribution. For "named" distributions, such as gaussian, there are often more efficient methods e.g. Box-Mueller. 
Inverse CDF computation can be expensive, because rarely you have the closed form expressions. There are, of course, approximations, such as this for Gaussian inverse CDF. In contrast, Box Muller approach uses certain properties of Gaussian distribution to reduce the algorithm to calculation of closed form expressions using square root, log and sine functions, all of which are implemented in the transcendental function modules of CPUs these days. The specialized random number generators will exploit special features of the distributions to get performance or precision edge over the generic method described above.
A very different approach is used in Markov Chain Monte Carlo methods. I'm not getting into this unless it's what you're looking for.
A: So here is an ECDF plot:

Here is the code used to make it:
library(latticeExtra)

y<- rnorm(n=1000,mean=0, sd=1)

ecdfplot(y)

Lets look at it.  It is comprised from 1000 samples of a standard normal distribution.  The x-axis is the magnitude of the sample, and the y-axis ranges from 0 to 1.  
The idea then is that you can make something that looks like this, but map from y to x.  Imagine that you drew 1000 samples from the uniform distribution, and then mapped them from y to x.  Your output would be a decent approximation of the standard normal distribution.
As long as you can make a continuous eCDF then you can use this method to approximately sample from the distribution.
So interpolation is a fair way to engage this.
p <- ecdf(y) #create function using data
fy <- p(sy)   #find values on sorted data (needed for interp)

yi <- sort(runif(n=1000)) #sample from uniform, sort

#use interpolation to map to distro of interest
x2 <- interp1(x=fy,y=sy,xi=yi,method="pchip")

#randomize if you need it
x2 <- sample(x2,length(x2),replace=FALSE)

#plot updated cdf and then compare
ecdfplot(x2)

So why Ecdf?  


*

*It is data driven.  If you have an analytic function for the
distribution then you can evaluate it at the points and make the
ecdf.  If you have only previous data and no analytic function then
you can make the ecdf.  If you want to hand-make your own
distribution using any valid mathematic method - as long as you end
up with numbers at the end, you can use this method to approximate
it.

*It is fast to code, and computers run very fast.  

*It works on data - and sometimes you can have cdfs for which analytic inverse does not exist.  


References:    


*

*http://blogs.sas.com/content/iml/2013/07/22/the-inverse-cdf-method.html

*https://en.wikipedia.org/wiki/Inverse_transform_sampling (aka Smirnov transform)

*http://www.mathworks.com/help/stats/examples/nonparametric-estimates-of-cumulative-distribution-functions-and-their-inverses.html
