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:


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
  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)

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)
        if(guesses[x] == chosenNum):
            numCorrect = numCorrect + 1
            closestNumber = min(guesses, key=lambda x:abs(x-chosenNum))
    averagePercentError = np.mean(percentErrors)
    print ("The average Percent Error is: " + str(round(averagePercentError,2)) + "%")

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

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")
        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
        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):
        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))

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.

  • 1
    \$\begingroup\$ What exactly is the basis for your methodology? Unless I'm missing something, this is just the german tank problem, for which there is a simple formula to estimate the total number of tanks from the max observed value. en.wikipedia.org/wiki/German_tank_problem#Example \$\endgroup\$ Dec 11, 2017 at 7:19
  • \$\begingroup\$ We are not supposed to know the answer to the problem and try to create a solution of our own. Also, our method has to use descriptive statistics not bayesian or frequentist inference. \$\endgroup\$ Dec 11, 2017 at 7:39
  • \$\begingroup\$ That makes perfect sense, thanks! I'd vote for maximum likelihood fitting myself, if possible: it would probably be quite a bit faster. I'm 95% sure though that it would be ruled out by those rules. \$\endgroup\$ Dec 11, 2017 at 12:59
  • \$\begingroup\$ It bugs me that you name all your functions starting with get even when their purpose is to produce output. \$\endgroup\$
    – Snowbody
    Dec 12, 2017 at 3:08

3 Answers 3


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)

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

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.

  • \$\begingroup\$ This was a great improvement! The program took 10.22 minutes for 10,000 guesses. \$\endgroup\$ Dec 11, 2017 at 1:11
  • \$\begingroup\$ There is still plenty of room for speedup. getNumGuessed is still quite slow for similar reasons. I'll update this answer once I have a good solution. \$\endgroup\$ Dec 11, 2017 at 1:13
  • \$\begingroup\$ Since STDarray is no longer a global variable, what do I put for STDarray.extend in the createListOfStandardDeviations(start,end)? \$\endgroup\$ Dec 11, 2017 at 1:33
  • \$\begingroup\$ At this point, the pace to gain speed is by in-lining createListOfStandardDeviationsand removing the loop around it through more vectorization \$\endgroup\$ Dec 11, 2017 at 1:36
  • \$\begingroup\$ I have an error saying that STDarray is not defined in the createListOfStandardDeviations(start,end). I do not know what to put there since we are defining STDarray from this function and I have no idea what we would be extending in that case. \$\endgroup\$ Dec 11, 2017 at 1:40

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 low1.

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)

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

Some simulation results

A closed form solution seems to be discussed here.

  • \$\begingroup\$ That's really annoying, as it makes vectorization way more difficult. Is there a way to make my version of createListOfStandardDeviations vectorized, fast, and correct? \$\endgroup\$ Dec 13, 2017 at 0:01
  • \$\begingroup\$ I don't see how my version makes vectorization any more or less difficult than the original version and how it should affect createListOfStandardDeviations. In both cases you get a 1-dimensional array tankSerialNumbers, only in my versions the numbers are guaranteed to be unique \$\endgroup\$
    – mkrieger1
    Dec 13, 2017 at 0:57
  • \$\begingroup\$ It is complicated because np.random.choice does not let you generate a 3d array where each element has a different maximum as I did in my version of createListOfStandardDeviations. This ability to vectorize led to much of the speed boost my solution provided. \$\endgroup\$ Dec 13, 2017 at 1:04
  • \$\begingroup\$ Then I seem to not quite understand what you're proposing. Anyway, the "offending" line tankSerialNumbers = np.random.randint(1, chosenNum + 1, size=numCapturedTanks) is outside the createListOfStandardDeviations function. \$\endgroup\$
    – mkrieger1
    Dec 13, 2017 at 1:12
  • \$\begingroup\$ How can I replace these two lines to pick without replacement maxes = np.broadcast_to(np.arange(start+1, end+1), (numCapturedTanks, 100, end - start)) tankSerialNumbersSimulated = np.floor(np.random.uniform(1, maxes)) \$\endgroup\$ Dec 13, 2017 at 1:13

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



(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)

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

  • \$\begingroup\$ Why are you randomizing the true number? I think you misunderstood the first OP. I'll clarify it in this post. \$\endgroup\$ Dec 11, 2017 at 2:31
  • \$\begingroup\$ I made it random so you could get more random simulations, that's why it's only initialized once, you can set it to whatever constant N that you want. \$\endgroup\$
    – mochi
    Dec 11, 2017 at 2:33
  • \$\begingroup\$ Ok. However, I don't get how you get one result from 10,000 guesses. I want 10,000 results from 10,000 guesses. \$\endgroup\$ Dec 11, 2017 at 2:46
  • \$\begingroup\$ Guesses is a list of num_guesses guesses, which in this case is a list of 10000 guesses. Where does it look like I'm only getting 1 result? Each iteration of the loop creates 1 guess from num_samples(in this case 100) simulated samples, by finding the max value of the simulated sample with the standard deviation closest to our observed standard deviation. \$\endgroup\$
    – mochi
    Dec 11, 2017 at 2:48
  • \$\begingroup\$ My bad with the 10,000 guesses thing. I had forgotten to write the calc_error function. \$\endgroup\$ Dec 11, 2017 at 2:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.