From 86128a6f00e11a06ea03e3b249d5f7db3af64a8a Mon Sep 17 00:00:00 2001 From: Holden Rohrer Date: Mon, 28 Mar 2022 19:13:51 -0400 Subject: more usable framework --- regression.py | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) create mode 100644 regression.py (limited to 'regression.py') diff --git a/regression.py b/regression.py new file mode 100644 index 0000000..b6ac7db --- /dev/null +++ b/regression.py @@ -0,0 +1,33 @@ +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) -- cgit