For-Loop with numpy array too slow
My toy example I build from your snippets:

import numpy as np


def throw_a_coin(N):
    return np.random.choice(['H','T'], size=N)


def make_throws(number_of_samples, sample_size):
    sample_probs = np.zeros(number_of_samples)
    for i in range(number_of_samples):
        num_heads = sum(throw_a_coin(sample_size) == "H")
        sample_probs[i] = float(num_heads) / sample_size
    return sample_probs

sample_sizes = [np.random.randint(1, 1e3) for idx in range(500)]

mean_of_sample_means = np.zeros(len(sample_sizes))
std_dev_of_sample_means = np.zeros(len(sample_sizes))

for i in range(len(sample_sizes)):
    prob = make_throws(200, sample_sizes[i])
    mean_of_sample_means[i] = np.mean(prob)
    std_dev_of_sample_means[i] = np.std(prob)  


cProfile (ordered by cumtime) reveals the problem (runs in 82.359 seconds):

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
 133/1    0.003    0.000   82.359   82.359 {built-in method builtins.exec}
     1    0.002    0.002   82.359   82.359 fast_flip.py:1(&lt;module&gt;)
   500    0.940    0.002   82.200    0.164 fast_flip.py:8(make_throws)
100000   79.111    0.001   79.111    0.001 {built-in method builtins.sum}
100000    0.058    0.000    2.148    0.000 fast_flip.py:4(throw_a_coin)
100000    1.455    0.000    2.090    0.000 {method 'choice' of 'mtrand.RandomState' objects}
100000    0.207    0.000    0.635    0.000 fromnumeric.py:2456(prod)
100000    0.026    0.000    0.429    0.000 _methods.py:34(_prod)
101500    0.409    0.000    0.409    0.000 {method 'reduce' of 'numpy.ufunc' objects}
     6    0.000    0.000    0.176    0.029 __init__.py:1(&lt;module&gt;)


There is a steep gap after buildins.sum, i.e. you spend most time there. We can use np.sum instead (pushing it down to 3.457 seconds):

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
 133/1    0.003    0.000    3.457    3.457 {built-in method builtins.exec}
     1    0.002    0.002    3.457    3.457 fast_flip.py:1(&lt;module&gt;)
   500    0.905    0.002    3.307    0.007 fast_flip.py:8(make_throws)
100000    0.046    0.000    1.869    0.000 fast_flip.py:4(throw_a_coin)
100000    1.287    0.000    1.823    0.000 {method 'choice' of 'mtrand.RandomState' objects}
201500    0.702    0.000    0.702    0.000 {method 'reduce' of 'numpy.ufunc' objects}
100000    0.172    0.000    0.536    0.000 fromnumeric.py:2456(prod)
100000    0.136    0.000    0.532    0.000 fromnumeric.py:1778(sum)
100000    0.023    0.000    0.378    0.000 _methods.py:31(_sum)
100000    0.021    0.000    0.364    0.000 _methods.py:34(_prod)
     6    0.000    0.000    0.178    0.030 __init__.py:1(&lt;module&gt;)


further we can replace the strings "H" and "T" with a Boolean and stay in numpy for longer (down to 1.633 seconds):

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
 133/1    0.003    0.000    1.633    1.633 {built-in method builtins.exec}
     1    0.001    0.001    1.633    1.633 fast_flip.py:1(&lt;module&gt;)
   500    0.101    0.000    1.485    0.003 fast_flip.py:11(make_throws)
100000    0.169    0.000    0.893    0.000 fast_flip.py:4(throw_a_coin)
100000    0.724    0.000    0.724    0.000 {method 'uniform' of 'mtrand.RandomState' objects}
100000    0.122    0.000    0.491    0.000 fromnumeric.py:1778(sum)
100000    0.024    0.000    0.354    0.000 _methods.py:31(_sum)
101500    0.334    0.000    0.334    0.000 {method 'reduce' of 'numpy.ufunc' objects}
     6    0.000    0.000    0.178    0.030 __init__.py:1(&lt;module&gt;)


next we can get rid of throw_a_coin and instead sample a numer_of_samples x sample_size array of uniformly distributed random numbers and threshold them. This also allows us to vectorize the for loop and stay in numpy even longer (0.786 seconds):

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
 133/1    0.003    0.000    0.786    0.786 {built-in method builtins.exec}
     1    0.003    0.003    0.786    0.786 fast_flip.py:1(&lt;module&gt;)
   500    0.053    0.000    0.634    0.001 fast_flip.py:4(make_throws)
   500    0.526    0.001    0.526    0.001 {method 'uniform' of 'mtrand.RandomState' objects}
     6    0.000    0.000    0.179    0.030 __init__.py:1(&lt;module&gt;)


Here is the code:

import numpy as np


def make_throws(number_of_samples, sample_size):
    # True == Heads
    throws = np.random.uniform(size=(number_of_samples, sample_size)) &gt; 0.5
    sample_probs = np.sum(throws, axis=1) / sample_size
    return sample_probs

sample_sizes = [np.random.randint(1, 1e3) for idx in range(500)]

mean_of_sample_means = np.zeros(len(sample_sizes))
std_dev_of_sample_means = np.zeros(len(sample_sizes))

for i in range(len(sample_sizes)):
    prob = make_throws(200, sample_sizes[i])
    mean_of_sample_means[i] = np.mean(prob)
    std_dev_of_sample_means[i] = np.std(prob)


At this point we could start to tackle the for loop, but it would start to get very micro optimized. One might consider aggregating the results in prob and then using np.mean(prob, axis=1) and np.std(prob, axis=1) outside the loop, but that only nets like 20ms so it's more of a personal preference thing.