Estimating the number of tanks based on a sample of serial numbers 2.0

Question

There was some confusion with the code I had posted in my previous version of this question and there was some good advice from @Oscar Smith.

The explanation is the same; please give improvements for speed in this new version:

#!/usr/bin/python

import numpy as np
import time

chosenNum = 429
numRuns = 10000
numCapturedTanks = 7
numGuessN = []
guesses = []
percentErrors = []
STDarray = []
start_time = time.time()
STDtimes = []

def getAverageStdTime(timetaken): # gets the average time it took to calculate standard deviations
  STDtimes.append(timetaken)
  if (len(STDtimes) == numRuns):
    print ("Average List of Standard Devations Generation Time: " + str(round(np.mean(STDtimes),2)) + " seconds")

def createListOfStandardDeviations(start,end):
    for y in range(start,int(end)):
        tankSerialNumbersSimulated = np.random.randint(1, y + 1, size=numCapturedTanks) #from Oscar Smith
        simulatedSTD = np.std(tankSerialNumbersSimulated)
        STDarray.append(simulatedSTD)

def getAllGuesses():
    print ("Your guesses are: " + str(guesses))

def getAvgPercentError():
    numCorrect = 0
    closestNumber = 0
    for x in range(len(guesses) - 1):
        percentError = '%.2f' % round(((np.abs(guesses[x] - chosenNum))/float(chosenNum) * 100), 2)
        percentErrors.append(float(percentError))
        if(guesses[x] == chosenNum):
            numCorrect = numCorrect + 1
        else:
            closestNumber = min(guesses, key=lambda x:abs(x-chosenNum))
    averagePercentError = np.mean(percentErrors)
    print ("The average Percent Error is: " + str(round(averagePercentError,2)) + "%")
    getAccuracy(numCorrect,closestNumber)

def getAccuracy(amountCorrect,closestNumberToActual):
    if (amountCorrect > 0):
        print ("You got the correct number " + str(amountCorrect) + " out of " + str(len(guesses)) + " times.")
    else:
        print ("Your closest number was: " + str(closestNumberToActual))
    getmode(guesses)

def getmode(inplist):
    dictofcounts = {}
    listofcounts = []
    for i in inplist:
        countofi = inplist.count(i) # count items for each item in list
        listofcounts.append(countofi) # add counts to list
        dictofcounts[i]=countofi # add counts and item in dict to get later
    maxcount = max(listofcounts) # get max count of items
    if maxcount ==1:
        print ("There is no mode for this dataset, values occur only once")
    else:
        modelist = [] # if more than one mode, add to list to print out
        for key, item in dictofcounts.items():
            if item ==maxcount: # get item from original list with most counts
                modelist.append(str(key))
        print ("Most guessed number(s):",' and '.join(modelist))
        return modelist

def getNumGuessed(givenSTD,maxNumber):
    minStd = min(STDarray, key=lambda x:abs(x-givenSTD)) #finds closest standard deviation to the given standard deviation
    for (z,this_std) in enumerate(STDarray):
        if(minStd == this_std): #find closest number to original standard deviation
            numGuessed = z + maxNumber
            return numGuessed

def main():
    print ("reached main")
    for runsRan in range(numRuns):
        tankSerialNumbers = np.random.randint(1, chosenNum + 1, size=numCapturedTanks) #from Oscar Smith
        NumSTD = np.std(tankSerialNumbers)
        highestTankSerial = np.mean(tankSerialNumbers) + 3*NumSTD
        maxNum = np.amax(tankSerialNumbers)
        print ("Tank Serial Numbers Generated")
        print ("Standard Deviation and Range Calculated")
        ListOfStandardDeviationsStartTime = time.time()
        for _ in range(100):
            del STDarray[:]
            if (maxNum - highestTankSerial < 0):
                createListOfStandardDeviations(maxNum,highestTankSerial)
            else:
                createListOfStandardDeviations(highestTankSerial,maxNum)
            numGuessN.append(getNumGuessed(NumSTD,maxNum))
        print ("Initial List of Standard Deviations Generated")
        print ("List of Standard Devations Generation took " + str(round(time.time() - ListOfStandardDeviationsStartTime,2)) + " seconds")

        guess = int(np.mean(numGuessN))
        print ("Guess Generated " + str(runsRan + 1))
        getAverageStdTime(float(time.time() - ListOfStandardDeviationsStartTime))
        guesses.append(guess)
    getAllGuesses()
    getAvgPercentError()

main()
print ("My program took " + str(round((time.time() - start_time)/float(60),2)) + " minutes to run")

Currently, the runtime is approximately 7.26 minutes for 1,000 runs. I want to get it to run 10,000 times and at this rate, it would take too long.

If anyone is confused by what the purpose is or any part of the code, please mention specifically what is confusing and I'll explain it.

To clarify: We are given 7 serial numbers based on a value n which is unknown (this number is given only to check how good our process of determining n is and to generate the 7 random serial numbers from it).

My process first finds the standard deviation of the given serial numbers. Then I find a limit that n definitely cannot exceed (three standard deviations above the mean) and a max from the given list as what it can't be below. I then simulate what random serial numbers would be generated from the predicted n from the range I found out. I take the standard deviation of each simulation and find which one is the closest to the standard deviation of the given serial numbers and store the corresponding guessed n. I do this x times (the more the better - I used 100) to get x guessed n 's. I take the mean of those guesses to get my final guess. I then find the percent error of my guess based on the actual number.

Answer 1

Some profiling at this point shows that almost all of the time is spent in createListOfStandardDeviations, so this is the next place to look for performance improvements. The way we can speed this up is by removing the for loop and making numpy do the work for us. In the old code, we calculated one simulatedSTD at a time, but it is possible to use a 2darray and compute the std of each row at the same time. Here is the updated code that does this.

def createListOfStandardDeviations(start,end):
    end = int(end)
    sim_num = end - start
    maxes = np.broadcast_to(np.arange(start+1, end+1), (numCapturedTanks, sim_num))
    tankSerialNumbersSimulated = np.floor(np.random.uniform(1, maxes))
    simulatedSTDs = np.std(tankSerialNumbersSimulated, axis=0)
    STDarray.extend(simulatedSTDs)

the main piece that is likely confusing here is the use of np.floor(np.random.uniform instead of np.randint. These do the same thing, and the switch is necessary because randint takes an int as the max, while uniform allows an ndarray instead. This change gives me a 35x speedup from before.

The next speedup we get is from making STDarray into an ndarray, and while we're at it, making it no longer a global variable. To do this, we make the last line of createListOfStandardDeviations be return np.std(tankSerialNumbersSimulated, axis=0), store it in main, and pass it into getNumGuessed. This change lets getNumGuessed become simply return maxNumber + np.argmin(np.abs(STDarray-givenSTD), axis=1) which is simpler and faster. This is a smaller boost than before, but it still gives about a 2x speedup.

At this point createListOfStandardDeviations is the slow point again, and the solution is yet more vectorization. At this point we are running this function in a loop 100 times, but it is much faster to let numpy do this. All we have to do is change the function to

def createListOfStandardDeviations(start,end):
    maxes = np.broadcast_to(np.arange(start+1, end+1), (numCapturedTanks, 100, end - start))
    tankSerialNumbersSimulated = np.floor(np.random.uniform(1, maxes))
    return np.std(tankSerialNumbersSimulated, axis=0)

and it will just work. It is important that getNumGuessed has the axis=0 in the argmin call or else it will not work properly. This gives another 2x speedup. Here is the main function after removing the loop.

def main():
    for runsRan in range(numRuns):
        tankSerialNumbers = np.random.randint(1, chosenNum + 1, size=numCapturedTanks) #generates initial tank serial numbers
        NumSTD = np.std(tankSerialNumbers)
        highestTankSerial = np.mean(tankSerialNumbers) + 3*NumSTD
        maxNum = np.amax(tankSerialNumbers)
        ListOfStandardDeviationsStartTime = time.time()
        start, end = min(maxNum, highestTankSerial), int(max(maxNum, highestTankSerial))

        STDarray = createListOfStandardDeviations(start, end)
        numGuessN.extend(getNumGuessed(STDarray, NumSTD, maxNum))

        guess = int(np.mean(numGuessN)) #store actual guess
        print ("Guess Generated " + str(runsRan + 1))
        guesses.append(guess) #add guesses to a list
    getAllGuesses()
    getAvgPercentError()

The last speedup worth bothering with is making numGuessN an ndarray. This will make finding it's mean much faster, and is easy enough to do. The three changes needed are changing the declaration to numGuessN = np.array([], dtype=np.int64), and changing the `extend call to

global numGuessN
numGuessN = np.concatenate((numGuessN, getNumGuessed(STDarray, NumSTD, maxNum)))

This yields another 2x speedup for me, At this point, we have gotten well over 200x over the original, and almost all of the time is spent either calculating std, mean or generating random numbers. It is possible that another factor of 2 could be gained from here, but it would be a pain, and it would probably require more effort than it's worth.

Answer 2

Drawing the serial numbers

tankSerialNumbers = np.random.randint(1, chosenNum + 1, size=numCapturedTanks)

This is the wrong way to generate the serial numbers, because it is possible that the same serial number is drawn more than once (this is random sampling with replacement):

>>> np.random.randint(1, chosenNum + 1, size=numCapturedTanks)
array([ 50,  10,  71, 244, 394, 375,  10]

You might not notice it, because the pool of possible serial numbers is much larger (chosenNum = 429) than the size of the sample you draw (numCapturedTanks = 7) and therefore the probability to draw the same number twice is low¹.

But the serial numbers really should be unique; you need random sampling without replacement.

This is done in plain Python using random.sample,

random.sample(range(1, chosenNum + 1), numCapturedTanks)

or in NumPy using np.random.choice with replace=False:

np.random.choice(range(1, chosenNum + 1), numCapturedTanks, replace=False)

This could be further simplified if you chose to use zero-based serial numbers (i.e. for 100 tanks, the serial numbers would be 0, …, 99 instead of 1, …, 100):

random.sample(range(chosenNum), numCapturedTanks)

np.random.choice(chosenNum, numCapturedTanks, replace=False)

---

¹Approximately 5% according to a quick Monte Carlo simulation I've run. (Source code)

A closed form solution seems to be discussed here.

Answer 3

You could vectorize your entire process in numpy.

From your previous post it looks like, from a given n number of tanks, you want to estimate a true N number of total tanks. The method is:

1. Generate m samples of n tanks
2. Calculate the std of these m samples
3. Find the sample with the closest std to your observed n tanks and use that as a guess to the true max, N
4. Repeat for a certain # of times to improve accuracy

So first, let's set some constants:

import numpy as np 

true_n = 429 # Can set whatever constant here
num_tanks = 7
num_guesses = 10000
num_samples = 100

Then a function for generating data:

def simulate_data(lower,upper,size):
    return np.array([np.random.randint(1,x+1,(num_tanks,size)) for x in range(lower, upper)])

Now we can generate our n sample and calculate relevant statistics:

obs = simulate_data(1,true_n,1)
obs_mean = np.mean(obs)
obs_std = np.std(obs)
obs_max = np.max(obs)
upper = int(obs_mean + 3*obs_std)

Finally we're ready to simulate data and run your algorithm:

min_guess, max_guess = np.sort([upper,obs_max])
guesses = []
for _ in range(num_guesses):  
    obs = simulate_data(1,true_n,1)
    obs_mean = np.mean(obs)
    obs_std = np.std(obs)
    obs_max = np.max(obs)
    upper = int(obs_mean + 3*obs_std)
    sim_data = simulate_data(min_guess, max_guess, num_samples)

    sim_std = np.std(sim_data,axis=1)
    min_std_offset = np.argmin(np.abs(sim_std - obs_std),axis=1)
    guess = min_guess + min_std_offset

    guesses.append(int(np.mean(guess)))

print("actual:",true_n)
print("guess:",np.mean(guesses))

(1) This combines your data generation with the total # of samples you want to make to generate all num_samples * num_tanks samples at once, in this case, 100 rows of 7 tank samples.

(2) Then calculates the standard deviation by row(axis=1) for all 100 samples.

(3) Next find the value that generated samples with the smallest std, this happens to be the starting value + whatever offset(index) the smallest value is at. Random samples are generated from (1,guess) where guess ranges from upper bound to the max value, or from max value to upper, whichever is larger.

(4) Stores the guessed value as a guess and repeats for 10000x.

This completes 10000 guesses for me in

25.7 s ± 1.01 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

If you had the memory for it, you could even add an additional dimension for num_guesses and skip the for loop.

You can also vectorize the way you calculate error in the same way:

def calc_error(guesses):
    diff = np.abs(guesses - true_n)
    perc_error = diff/true_n*100
    num_correct = sum(guesses == true_n)
    closest_guess = guesses[np.argmin(diff)]

    print("num correct:",num_correct)
    print("avg percent error: %.2f" % np.mean(perc_error))
    print("closest guess was:",closest_guess)

Then you just need to convert your guesses list into a numpy array and call the error function on it:

guesses = np.array(guesses)
calc_error(guesses)

This follows your procedure now I believe, and from what I can tell finishes in a fraction of the time of your current solution.