Simulation of Central Limiting Theorem (CLT)

Question

The following is the simulation of Central Limiting Theorem from Statistics:

class MyRandom:
 m = 4294967296;# modulus
 a = 1664525; # multiplier
 c = 1013904223; # increment
 seed = 1; #initial seed
 def NextInt(self):
 self.seed = (((self.a * self.seed + self.c) % self.m));
 return self.seed;
 def NextIntMm(self, mins, maxs):
 temp = self.NextInt();
 ddd = float(temp) / float(self.m);
 returns = int((maxs - mins) * ddd + mins);
 return returns;
 def NextDouble(self):
 temp = self.NextInt();
 return float(temp) / float(self.m);
 def NextDoubleMm(self, mins, maxs):
 temp = self.NextInt();
 fraction = float(temp) / float(self.m);
 return (maxs - mins) * fraction + mins;
#################################################################
N = 999999; 
r = MyRandom();
population = [];
population_sum = 0;
for i in range(0, N):
 x = r.NextDouble();
 population.append(x);
 population_sum = population_sum + x;
population_mean = population_sum / len(population);
print("Population mean = ", population_mean);
##################################################################
rnd_num_count = 30;
sample_size = 12;
sample_means = [];
sum_of_sample_mean = 0;
for i in range(0, rnd_num_count):
 sums = 0;
 for j in range(0, sample_size):
 index = r.NextIntMm(0, N - 1);
 element = population[index];
 sums = sums + element;
sample_mean = sums / sample_size;
sum_of_sample_mean = sum_of_sample_mean + sample_mean;
sample_means.append(sample_mean);
print("Mean of sample mean = ", sum_of_sample_mean/ len(sample_means));

Here, Normal Distribution of population is used.

Can you review this and point out any possible improvement needed?

N.B. Defining a random number generator is part of the requirement.

• Graipher: You are aware of the random module in the standard library, right? If so, consider adding the reinventing-the-wheel tag.
• Graipher: I guess you have an indentation error in your second example, the second-to-last block should probably under the outer for loop.
• user366312: @Graipher, No, my program is running smoothly in Spyder.
• Graipher: That may be the case, but it will probably not do the right thing because sample_means would contain only the last sample mean and sums will be reset every loop iteration without being used, so sample_mean will only use the sums from the last iteration.

Answer 1

Even though you are implementing your own random number generator, I would stick to the interface defined in the standard library module random. You could have a random, a uniform, a randint (or randrange if it is inclusive on the end) and a choices method.

Your class should probably also allow setting the seed on creation.

Python has an official style-guide, PEP8, which recommends not using trailing ; and using lower_case both for variables as well as functions/methods.

Even in Python 2 (which you should no longer use), if either the numerator or the divisor is a float, the result is also a float with the division operator /, so you only need to cast one of them. In Python 3 / is always float division and // is always integer division, so no need for the casts to float.

You should always add a dosctring to your code explaining how to use it.

For range, the start is implicitly 0. Also _ is conventionally used for unused loop variables.

Wrap the code that uses your class in a if __name__ == "__main__": guard to allow importing from this script without it running.

class Random:
 """A pseudo-random number generator using a LCG."""
 m = 4294967296 # modulus
 a = 1664525 # multiplier
 c = 1013904223 # increment
 def __init__(self, seed=1):
 """Initialize the pseudo-random number generator with a seed (default: 1)."""
 self.state = seed
 def randint(self, a=None, b=None):
 """Return either a random integer between `0` and `self.m`
 or between `a` and `b` (exclusive?)."""
 if a is None and b is None:
 self.state = (self.a * self.state + self.c) % self.m
 return self.state
 elif a is not None:
 if b is not None:
 return int(self.uniform(a, b))
 else:
 return int(self.uniform(0, a))
 else:
 raise ValueError("Need to also set `a` if you set `b`.")
 def random(self):
 """Return a random float between 0 and 1 (exclusive?)."""
 return float(self.randint()) / self.m
 def uniform(self, a, b):
 """Return a random float between `a` and `b` (exclusive?)."""
 return (b - a) * self.random() + a
 def choices(self, population, k=1):
 """Return a k sized list of population elements chosen with replacement."""
 n = len(population)
 return [population[self.randint(n - 1)] for _ in range(k)]

The usage of this is similar to your code, just with renamed methods.

You can also greatly shorten your examples using list comprehensions and a utility function mean:

random = Random()
def mean(x):
 return float(sum(x)) / len(x)
if __name__ == "__main__":
 N = 999999 
 population = [random.random() for _ in range(N)]
 print("Population mean = ", mean(population))
 rnd_num_count = 30
 sample_size = 12
 sample_means = [mean(random.choices(population, k=sample_size))
 for _ in range(rnd_num_count)]
 print("Mean of sample mean = ", mean(sample_means))

This script has the advantage that you can easily test it with the standard library random module by replacing the line random = Random() with import random. This way you can see if your class produces random numbers somewhat correctly.

---

If implementing so many different random functions gets tedious, you can also subclass random.Random:

Class Random can also be subclassed if you want to use a different basic generator of your own devising: in that case, override the following methods: random(), seed(), getstate(), and setstate(). Optionally, implement a getrandbits() method so that randrange() can cover arbitrarily large ranges.

This of course means that you are using the standard library and not your own implementations as much. In that case you probably also want to use statistics.mean.