# Suppose you fund 100,000 scientific experiments. You go off, sample
# the impact of 100 of them. How does the impact of the total compare
# to the impact of the sample?
#
# If you're like most people, you'd guess it's about 1,000 times as
# large. This is often true. But it's also often false: when the
# sample is tail dominated. The total set will contain many more
# extreme outliers, which will tend to greatly increase the average
# impact.
#
# This phenomenon makes it difficult to compare intervention and
# control strategies for funding.
#
# This quick-and-dirty script makes this comparison explicit, as
# described in the essay: Michael Nielsen and Kanjun Qiu, "A Vision of
# Metascience: An Engine of Improvement for the Social Processes of
# Science", https://scienceplusplus.org/metascience/index.html, San
# Francisco (2022).
import math
import numpy as np
import random
random.seed(314159) # for repeatability
# power law
def power_sample(alpha, delta):
return random.random()**(-1/(alpha-1))-delta
def repeated_power_sample(n, alpha, delta):
return [power_sample(alpha, delta) for j in range(n)]
def median(alpha, delta):
return 2**(1/(alpha-1))-delta
# lognormal
def repeated_lognormal_sample(n, mu, sd):
return [math.exp(random.gauss(mu, sd)) for j in range(n)]
n_trial = 100 # size of the sample
n = 100000 # total size
alpha = 2.3 # shape parameter for the power law for the intervention
delta = 1 # displace intervention power law value by this amount,
# so minimum value is 0 (not 1)
mu = 0 # mean for lognormal for the control
sd = 1.4 # sd for lognormal for the control
num_comparisons = 1000 # We'll make our comparison many times
def comparison(verbose):
I = [y for y in repeated_power_sample(n, alpha, delta)]
C = repeated_lognormal_sample(n, mu, sd)
I_trial = sum(I[:n_trial])
C_trial = sum(C[:n_trial])
I_t = sum(I)
C_t = sum(C)
if verbose:
print("The control distribution is a lognormal C with mu = {}, sd = {}, and median {}".format(mu, sd, math.exp(mu)))
print("The intervention is a power law with shape alpha = {}, and displaced by an amount {},\nso the range is 0...infinity. The median is {}".format(alpha, delta, median(alpha, delta)))
if median(alpha, delta) < math.exp(mu):
print("The intervention I has a smaller median than the control C")
else:
print("The intervention I has a large median than the control C")
print("I trial has average impact {}".format(I_trial/n_trial))
print("C has average impact {}".format(C_t/n))
if I_trial/n_trial < C_t/n:
print("The I trial has lower average impact than C overall")
else:
print("The I trial has higher average impact than C overall")
print("I done at scale has average impact {}".format(I_t/n))
if I_t/n < C_t/n:
print("I done at scale has lower average impact than C overall")
else:
print("I done at scale has higher average impact than C overall")
return (median(alpha, delta) < math.exp(mu)) and \
(I_trial/n_trial < C_t/n) and (I_t/n > C_t/n), I_trial/n_trial, I_t/n
verbose = False # no commentary as we go, but useful for building insight to turn on
comparisons = [comparison(verbose) for j in range(num_comparisons)]
occ = sum(int(x[0]) for x in comparisons)
print("{} of {} comparisons have the claimed property".format(occ, num_comparisons))
shifts = [x[2]/x[1] for x in comparisons]
average_shift = sum(shifts)/num_comparisons
print("The average shift is {}".format(average_shift))