Hi all, this question doesn't really have to do with C, but it is related to simulation, and without solving this first, I can't do the random variate generator.

I'm supposed to generate random variates for p(x) ~ pi - x + x^2 + sin(x), for 0 <= x <= pi.

I got the normalization factor as 6 / (12 + 3*pi^2 + 2*pi^3), which gives me my f(x) as

f(x) = [6 / (12 + 3*pi^2 + 2*pi^3)] * (pi - x + x^2 + sin(x)).

I then run into problems decomposing the f(x). I decomposed it into the following equation:

f(x) = 6*pi / a - (3*pi^2 / a) * f_1(x) + (2*pi^3 / a) * f_2(x) + (12 / a) * f_3(x)
where a = (12 + 3*pi^2 + 2*pi^3), f_1(x) = 2x/pi^2, f_2(x) = 3x^2 / pi^3, f_3(x) = ½ sinx.

I reasoned that f_1 is proportional to x, f_2 to x^2, and f_3 to sin(x). I then normalized f_1, f_2 and f_3, obtaining the respective individual distributions.

However, I'm confused by the presence of the constant term 6*pi / a and the negative sign in front of f_1. What should I do with them?

I'm really confused with the whole composition method, and I would greatly appreciate it if someone could point me to a website that gives a clear explanation of how to decompose the distributions.

Thank you.

Regards,
Rayne

Well, what I suggest is don't decompose
[6 / (12 + 3*pi^2 + 2*pi^3)] * (pi - x + x^2 + sin(x))

The term [6 / (12 + 3*pi^2 + 2*pi^3)] is a constant. Evaluate it first and hold it as a constant to plug into the equation when needed.

__________________

Got a cough? Go home tonight and eat a whole box of Ex-Lax. Tomorrow, you'll be afraid to cough.
-- Pearl Williams

I'm sorry, but I still don't understand. I have to use composition method. If I don't decompose my normalized p(x), how do I generate random variates?

Hi all,
I'm supposed to generate random variates using the M(RT)^2 method, for p(x) ~ pi - x + x^2 + sin(x), for 0 <= x <= pi. I normalised p(x) to get p(x) = [6/(12 + 3pi^2 + 2pi^3)]*(pi - x + x^2 + sin(x)).

I wrote the code below, to find out which value my parameter "ar" should be in order for the acceptance ratio to be 2/3. However, I kept getting an acceptance ratio of as high as 99.5%, regardless of my "ar" value. Also, I noticed that the random variates x_old and x_new just kept on either increasing or decreasing till they are very far off the range of [0, pi]. My probabilities of x_old and x_new (p_xold and p_xnew respectively) also went into the hundreds range the more points I tried to take. This is obviously wrong, isn't it?

What is wrong with my code?

Plesae advise.


Thank you.

Regards,
Rayne

Maybe I'm missing something, because I don't see how your technique is supposed to work. In particular, I don't think what you are attempting implements what is usually called the "composition method"

I'll try to explain what I think is meant by the "composition method":

As I see it, if you want to generate variates from a given probability density function (pdf), one way is by inversion. This requires that you get samples from a uniform pdf and plug them into the inverse function of f(x). That's not easy with your the given pdf, so you consider the composition method:

If f(x) is not easily invertable, then you can try to express f(x) as a sum of terms, where each term (with suitable scale factor) is a pdf, and the regions covered by the terms don't overlap. If you can do this, and you can find the inverse function of each the terms, then there is a way to get several samples (depending on the number of terms) from a uniform pdf and create a sample from f(x). I won't go into the exact calculations required to generate the variates after you come up with the terms that this method needs.

The trick, of course, is picking the terms that meet these qualifications:

1. The sum of the terms is equal to the given pdf.

2. Each term is a constant scale factor times a pdf function in some form

3. You can find the inverse function of each term

4. The regions defined by each term do not overlap with regions from the other terms.

If you can do this, then there is a way to generate variates from the given pdf by combining several independent samples obtained from a uniform distribution.

I don't see how your f(x) can be broken up into terms like this. Note that, in general, if you have f(x) = p1 * g1(x) + p2 *g2(x), the inverse of f(x) is not necessarily equal to the sum of the inverses of the terms on the right hand side (unless the regions covered by the terms don't overlap).

So, I don't see where your f(x) is a suitable candidate for this approach, since I don't see how to get terms that satisfy the conditions that I stated. Note that just because I say that, "I don't see...," doesn't mean that it can't be done. Maybe someone else can show us how.

I suggest that if your pdf doesn't lend itself to methods that involve inverting functions, that you consider the Acceptance-Rejection method.

Regards,

Dave

I think you were replying to my second problem...
OK, I'll try to clarify the situation. I was given the same p(x) ~ pi - x + x^2 + sin(x) for 0 <= x <= pi. I then had to generate random variates using the composition method, acceptance-rejection method, AND the M(RT)^2 (Metropolis Algorithm) method.

I had already completed the acceptance-rejection method, but am still stuck at the composition method and the M(RT)^2 method. My first post was about the composition method, and my second post (the one with the inserted code) was about the M(RT)^2 method.

I'm stuck at the composition method because I had trouble decomposing the normalized p(x). I am stuck at the M(RT)^2 method because I couldn't find the correct value of the parameter "ar" for generating random new x-points so as to get an acceptance rate of 67%. My acceptance rate is always above 95% for some reason.

So basically I need help in using these 2 methods to generate random variates for the p(x).

Thank you.

Regards,
Rayne

You are correct: I read your first post (about the composition method) and looked at the code for your second post. Since I can't see how to help with the composition, and you already have solved the acceptance-rejection method, and I have never actually used the M(RT)^2 methid, I don't have any further suggestions. I am sorry that I can't be more helpful.

Regards,

Dave

That's all right, thank you so much for taking the time to help me anyway.

Regards,
Rayne