Implement density function

Question

I am going through my book , it states that " Write a sampling algorithm for this density function"

y=x^2+(2/3)*x+1/3; 0 < 𝑥 < 1

Or I can use Monte Carlo? Any help would be appreciated!

Answer 1

I'm assuming you mean you want to generate random x values that have the distribution specified by density y(x).

It's often desirable to derive the cumulative distribution function by integrating the density, and use inverse transform sampling to generate x values. In your case the the CDF is a third order polynomial which doesn't factor to yield a simple cube-root solution, so you would have to use a numerical solver to find the inverse. Time to consider alternatives.

Another option is to use the acceptance/rejection method. After checking the derivative, it's clear that your density is convex, so it's easy to create a bounding function b(x) by drawing a straight line from f(0) to f(1). This yields b(x) = 1/3 + 5x/3. This bounding function has area 7/6, while your f(x) has an area of 1, since it is a valid density. Consequently, 6/7 of points generated uniformly under b(x) will also fall under f(x), and only 1 out of 7 attempts will fail in the rejection scheme. Here's a plot of f(x) and b(x):

Since b(x) is linear, it is easy to generate x values using it as a distribution after scaling by 6/7 to make it a valid distribution function. The algorithm, expressed in pseudocode, then becomes:

function generate():
  while TRUE:
    x <- (sqrt(1 + 35 * U(0,1)) - 1) / 5     # inverse CDF transform of b(x)
    if U(0, b(x)) <= f(x):
      return x
  end while
end function

where U(a,b) means generate a value uniformly distributed between a and b, f(x) is your density, and b(x) is the bounding function described above.

I implemented the algorithm described above to generate 100,000 candidate values, of which 14,199 (~1/7) were rejected, as expected. The end results are presented in the following histogram, which you can compare to f(x) in the plot above.