Member Avatar for highflyer8

Hi, I don't understand the need for the while statement. Any help would be greatly appreciated.

double michel()
{
	double x, y;
	do                                // The loop computes points on the graph of y against x.
	{
		x = 53.0*rand()/RAND_MAX;     // x is the momentum. x is randomly generated. Min x is 0 MeV/c. Max x is 53 MeV/c.
		y = x/53.0;                   // y is the probability distribution function. Min y is 0.       Max y is 1.
	}                                 // y is a triangular distribution. Hence, y = (1/53)*x.
	while (1.0*rand()/RAND_MAX > y);  // Because we want to start with the WHYYYYYYYYYYY????????????
	return x;
}
  • 3 Contributors
  • 2 Replies
  • 143 Views
  • 9 Hours Discussion Span
  • Latest Post Latest Post by template<>

Dani AI

Generated

A do/while is appropriate here because the loop body must produce at least one candidate before any acceptance test runs. The pattern in the thread is an accept–reject sampler: the body generates a candidate x and computes an acceptance probability (proportional to x), and the loop condition draws a fresh uniform random number and repeats while that draw exceeds the acceptance probability. That is why the check sits after the body rather than before it.

’s suggested comparison (comparing x/53.0 to y) misunderstands the intent: in the original code y is derived from x, so comparing them is tautological and will not implement acceptance. ’s nudge to explain the goal is on point — the aim is to produce values whose density is proportional to x on [0, 53] (a simple triangular-like distribution).

A simpler, exact approach is to use the inverse CDF instead of rejection sampling. For density proportional to x on [0,L] the CDF is (x/L)^2, so x = L * sqrt(U) for U ~ Uniform(0,1). Example using modern C++ random facilities:

#include <random>
#include <cmath>

std::mt19937 rng{std::random_device{}()};
std::uniform_real_distribution<double> u(0.0, 1.0);
double x = 53.0 * std::sqrt(u(rng));

Practical tips: prefer <random> over rand(), seed your engine intentionally, remember the do/while semicolon is required, and note that rejection here would on average take about 2 tries (mean acceptance ~0.5).

Recommended Answers

All 2 Replies

Member Avatar for Saith

Neither do I, really. I don't see how a random number should be existing in this loop, and then an X value to be returned. I can understand if you want to compare X > Y and change the while expression to

while (x / 53.0 > y);

so that a more valid x is returned when compared to a y value.


(With the original while loop). Imagine if the first loop, x was 50 and the random expression in the while (1 > .6), the loop would go again. Say by happenstance, x was 50 (again) and the random expression in the while loop was now (.6 > 1). The statement would turn false and now you would return x with value 50. To me this seems like a logical error. Granted, I don't know much more of what is going on other than the code that was given.

If the while loop was changed to

while (x / 53.0 > y);

x = 50.0 and then in while statement you'd get (50.0/53 > y), I can then understand why you would return the x value from the expression of the y value.

You may also want to check out this for Probability Distribution on wiki.

http://en.wikipedia.org/wiki/Cumulative_distribution_function


Finally, I could easily be wrong right from the get-go. I have never taken statistics but would be interested in seeing the answer.

Edited by Saith because: n/a
Member Avatar for template<>

The code does not make any sense. Would be helpful to describe what you are trying to do in words. For example what are you trying to do with the while loop? What is the problem you're trying to solve?

Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, networking, learning, and sharing knowledge.