Python code translation issue or language difference in random number generation

Question

I'm new to Python and have no experience with QBasic. I was running a simulation in Python that came up with theoretically wrong values. I then ran it in QBasic and came up with theoretically predicted values.

Here are the test cases. I'm only counting the probability P(0.9<%Y<=1.8) so the count has to fall within those values. The 1-random.random() was only for that case, when I tried using that for all the cases they still came up with the wrong values. Here's the theoretical outcomes and you can see how it's different:

y~u(0,1) = 0.575

y~exp(2) = 0.3371

x1~u(0,1) x2~u(0,2) = 0.4475

P(y=0.25)=0.8 P(y=1.5)=0.2 = 0.32

---

In Python the simulation code is:

def test():
  x1,x2,c = 0.0,0.0,0.0
  for i in range(10000):
    if random.random()< 0.8:
      x1 += 0.25
    else:
      x2 += 1.5
    y = x1 + x2
    if y>0.9 and y<=1.8:
      c = c + 1
  return x1,x2,c

print "test: ",test()

def sim(a,b):
  #pyab1 = sum([a for a in a if a>0.9 and a<=1.8])/10000
  #pyab2 = sum([b for b in b if b>0.9 and b<=1.8])/10000
  #print "*****",float(pyab1+pyab2)
  #print a+b
  #array1 = [[a],[b]]
  array1 = a+b
  #array1.extend(a)
  #array1.extend(b)
  #c = 0
  #for y in array1:
    #if y>0.9 and y<=1.8:
      #c = c + 1
  pyab = sum([y for y in array1 if y>0.9 and y<=1.8])/10000
  print("P(a < x <= b) : {0:8.4f}".format(pyab))

Here's the Python output followed by the values it's supposed to give, but this shows how far off the results are.

case 1: P(a < x <= b) : 0.7169 #should be 0.575 
case 2: P(a < x <= b) : 0.4282 #should be 0.3371 
case 3: P(a < x <= b) : 0.5966 #should be 0.4475 
case 4: P(a < x <= b) : 0.5595 #should be 0.32 

---

In QBasic the simulation code is:

Case 1:

RANDOMIZE
FOR i = 1 TO 10000
    X1 = RND(1)
    X2 = RND(1)
    Y = X1+X2
    IF (Y>0.9) AND (Y<=1.8) THEN C=C+1
NEXT i

PRINT C/10000

Case 2:

RANDOMIZE
FOR i = 1 TO 10000
    X1 = (-0.5)*(LOG(1-RND(1)))
    X2 = (-0.5)*(LOG(1-RND(1)))
    Y = X1+X2
    IF (Y>0.9) AND (Y<=1.8) THEN C=C+1
NEXT i

PRINT C/10000

Case 3:

RANDOMIZE
FOR i = 1 TO 10000
    X1 = RND(1)
    X2 = RND(1)*2
    Y = X1+X2
    IF (Y>0.9) AND (Y<=1.8) THEN C=C+1
NEXT i

PRINT C/10000

Case 4:

RANDOMIZE
FOR i = 1 TO 10000
    X14 = RND(1)
    X24 = RND(1)
    IF (X14<0.8) THEN X41=0.25 ELSE X41=1.5
    IF (X24<0.8) THEN X42=0.25 ELSE X42=1.5
    Y = X1+X2
    IF (Y>0.9) AND (Y<=1.8) THEN C=C+1
NEXT i

PRINT C/10000

Here's the QBasic output, which shows how this is actually getting the right results.

case 1: P(a < x <= b) :  0.5715 
case 2: P(a < x <= b) :  0.3371 
case 3: P(a < x <= b) :  0.4413 
case 4: P(a < x <= b) :  0.3213 

All of the above code works for me without error. I don't see any differences in the algorithm used to get the values. Not sure if Python generates numbers differently from QBasic, and if that accounts for the reason behind this behavior.

I'm new to both of these languages but QBasic seems very primitive, and it would seem more likely that Python would get the right answer and QBasic the wrong one. But the opposite is happening. It doesn't appear to be something related to any difference in code. On translation, they both seem to be saying the same thing.

I'm interested in the reason why they give the 2 different outcomes. I'm more interested in why Python is giving the wrong answers and QBasic is giving the right ones.

Answer 1

Your Python code is totally wrong. What I think you want it to do is the following:

• Take two arrays, a and b, each containing 10000 numbers generated by some random function. (Equivalently, each an array of 10000 samples from data following some given distribution.)
• Pair up the values into 10000 pairs, with each pair taking an element from a and an element from b.
• Take the sum of each pair.
• Count how many of those 10000 pair-sums lie between 0.9 and 1.8
• Divide the above count by 10000 to get the probability that any given pair of data drawn from these distributions will sum to between 0.9 and 1.8, and print that probability.

However, your sim(a, b) function is doing something wildly different. Basically, what you're actually doing is:

• Concatenate the two 10000 element arrays, forming a 20000-element array of their elements.
• Take the sum of any elements in this new 20000-element array that are greater than 0.9.
• Divide that sum by 10000 and print it.

This algorithm bears no resemblance to anything in your Q-Basic code.

If I've understood your problem properly, I think what you want your sim function to be is this:

def sim(x_sample, y_sample):

    count = 0
    for i in range(10000):
        if 0.9 <= x_sample[i] + y_sample[i] <= 1.8:
            count += 1

    probability = count/10000.0

    print("P(a < x <= b) : {0:8.4f}".format(probability))

(There are almost certainly more elegant and Pythonic ways to implement the above function, but this way should be easy for a Python newbie to understand.)

Here are the results of tests I ran in the interpreter for your cases 1 to 3, as I've understood them from the QBasic program. I haven't included a version of test 4 because I didn't understand the QBasic code for test 4. The results for the first three tests are what you said they should be, though.

>>> from random import random
>>> 
>>> sim([random() for i in range(10000)],
...     [random() for i in range(10000)])
P(a < x <= b) :   0.5746
>>>      
... from math import log
>>> 
>>> sim([-0.5*log(1-random()) for i in range(10000)],
...     [-0.5*log(1-random()) for i in range(10000)])
P(a < x <= b) :   0.3405
>>>     
... sim([random() for i in range(10000)],
...     [2*random() for i in range(10000)])
P(a < x <= b) :   0.4479