import numpy as np

# data = [(0,0), (1,1), (2, 4), (3, 9), (4, 16) ] ; t = 3
def velvar(data, t):
    n = len(data) # number of entries
    X = np.array([(datum[0]**2, datum[0], 1) for datum in data])
    Xt = np.transpose(X) # X^T
    Y = [datum[1] for datum in data]
    M = np.matmul( np.linalg.inv(np.matmul(Xt, X)), Xt ) # solution = My
    est = np.matmul(M, Y);
    P = np.matmul(X, M) # \hat y = Py
    R = np.subtract(np.identity(n), P) # residual = Ry
    bias = sum([sum([entry**2 for entry in row]) for row in R])
    Ry = np.matmul(R, Y);
    S2 = np.dot(Ry, Ry); # S^2 estimator
    variance = S2/bias
    v = 2*M[0]*t + M[1] # V = v\cdot y
    velocity = 2*est[0]*t+est[1]
    variancevel = np.dot(v,v)*variance
    return (np.abs(velocity), variancevel)

def combine_estimates(velvars):
    avg = np.average(list(a[0] for a in velvars))
    var = (np.sum(a[1] + a[0]**2 for a in velvars) - avg**2)/len(velvars)
    return((avg, var))

if __name__ == '__main__':
    tot = (velvar([(0,0), (1,1), (2, 4), (3, 9), (4, 17)], 3),
           velvar([(0, 0), (1,2), (2, 8), (3, 18), (4, 32)], 3))
    n = 5
    k = 2
    print(sum([datum[1]**2+n*datum[0]**2 for datum in tot])
        - n*k*sum([datum[0]/k for datum in tot])**2)