4822.sagews

def statistical_power(x, n1, n2, alpha):
    """
    Return a statistical power.
    
    INPUT:
        x --- (myu2 - myu1)/sigma
        n1, n2 --- sample sizes
        alpha --- rate of the type I error
    OUTPUT:
        A statistical power
    """
    
    T = RealDistribution("t", n1 + n2 - 2)
    F = T.cum_distribution_function
    c = sqrt((n1*n2)/(n1 + n2))
    t = -T.cum_distribution_function_inv(alpha/2)
    return F(c*x - t) + F(-c*x - t)

def fast_find_min(alpha=0.05, beta=0.05, delta=0.5):
    """
    Return a minimum sample size pair if two sample sizes are *equal*.
    INPUT:
        alpha --- ratio of the type I error
        beta --- ratio of the type II error
        delta --- permissible difference
    OUTPUT:
        A minimum sample size.
    """
    
    f = RealDistribution("gaussian", 1).cum_distribution_function_inv
    u_alpha = -f(alpha/2)
    u_beta = -f(beta)
    b = floor(2*(((u_alpha + u_beta)/delta)^2))
    while statistical_power(delta, b, b, alpha) < 1 - beta:
        b *= 2
        
    ## binary search
    a = 1
    while a < b:
        mid = int((a + b)/2)
        if statistical_power(delta, mid, mid, alpha) < 1 - beta:
            a = mid + 1
        else:
            b = mid
    return a

for delta in [1, 0.5, 0.1, 0.05, 0.01, 0.005, 0.001]:
    print delta
    time print fast_find_min(alpha=0.05, beta=0.05, delta=delta)

1 28 Time: CPU 0.01 s, Wall: 0.01 s 0.500000000000000 106 Time: CPU 0.02 s, Wall: 0.02 s 0.100000000000000 2602 Time: CPU 0.02 s, Wall: 0.02 s 0.0500000000000000 10398 Time: CPU 0.02 s, Wall: 0.02 s 0.0100000000000000 259896 Time: CPU 0.03 s, Wall: 0.03 s 0.00500000000000000 1039578 Time: CPU 0.04 s, Wall: 0.04 s 0.00100000000000000 25989420 Time: CPU 0.04 s, Wall: 0.04 s

alpha = 0.05
beta = 0.05
a = 0.5
b = 1
f = RealDistribution("gaussian", 1).cum_distribution_function_inv
u_alpha = -f(alpha/2)
u_beta = -f(beta)
u = 2*((u_alpha + u_beta)^2)
answers = points([ (delta, fast_find_min(alpha=alpha, beta=beta, delta=delta)) for delta in srange(a, b, a/10.) ], size=15)
approx = plot( u/(x^2), (x, a, b), color="red", thickness=2)
show(answers + approx)
(answers + approx).save("slight_difference.eps")

def laymans_find_min(cost, budget, epsilon, alpha=0.05, beta=0.05, delta=0.5):
    """
    Return a minimum sample size pair in a sense of cost function.
    NOTE:
        This may take some time;
        because it searchs MANY.
    INPUT:
        cost --- linear function N^2 -> R such that c1 n1 + c2 n2
        budget --- maximal cost
        epsilon --- region
        alpha --- ratio of the type I error
        beta --- ratio of the type II error
        delta --- permissible difference
    OUTPUT:
        A pair of natural numbers that makes cost function minimum.
    """
    
    f = RealDistribution("gaussian", 1).cum_distribution_function_inv
    u_alpha = -f(alpha/2)
    u_beta = -f(beta)
    c1 = cost((1, 0))
    c2 = cost((0, 1))
    N1 = int(float(budget)/c1)
    N2 = int(float(budget)/c2)
    u = ((u_alpha + u_beta)/delta)^2
    domain = [ (n1, n2) for n1 in range(1, 1 + N1) for n2 in range(1, 1 + N2)
        if n1 + n2 > 2 and
        c1*n1 + c2*n2 <= budget and
        u < (n1*n2)/(n1 + n2) < u + epsilon and
        statistical_power(delta, n1, n2, alpha) >= 1 - beta]
    if domain == []: return IOError, "Too strict condition. No answer found."
    else: return min(domain, key=cost)

alpha = 0.05
beta = 0.05
delta = 0.5
epsilon = 2
cost = lambda (n1, n2): n1 + 2*n2
budget = 500
f = RealDistribution("gaussian", 1).cum_distribution_function_inv
u_alpha = -f(alpha/2)
u_beta = -f(beta)
c1 = cost((1, 0))
c2 = cost((0, 1))
N1 = int(float(budget)/c1)
N2 = int(float(budget)/c2)
u = ((u_alpha + u_beta)/delta)^2
domain = points([ (n1, n2) for n1 in range(1, 1 + N1) for n2 in range(1, 1 + N2)
        if n1 + n2 > 2 and
        c1*n1 + c2*n2 <= budget and
        u < (n1*n2)/(n1 + n2) < u + epsilon and
        statistical_power(delta, n1, n2, alpha) >= 1 - beta])
show(domain)

for n in range(2, 10):
    time print n, laymans_find_min(lambda (n1, n2): n1/float(n) + (n - 1)*n2/float(n), 150, 2, delta=0.5)

2 (106, 106) Time: CPU 1.71 s, Wall: 1.71 s 3 (127, 91) Time: CPU 1.91 s, Wall: 1.91 s 4 (141, 85) Time: CPU 2.22 s, Wall: 2.22 s 5 (150, 82) Time: CPU 2.78 s, Wall: 2.78 s 6 (166, 78) Time: CPU 2.89 s, Wall: 2.89 s 7 (187, 74) Time: CPU 3.51 s, Wall: 3.51 s 8 (187, 74) Time: CPU 3.59 s, Wall: 3.59 s 9 (201, 72) Time: CPU 4.24 s, Wall: 4.24 s