Privacy Preserving Recursive Least Squares

17 days ago by ktjell12@student.aau.dk

##################################################################################################### # This cell implementes most of the SMPC funtions, used in the paper: # # Privacy Preserving Recursive Least Squares Solutions by Katrine T., Rafael W., Ignacio C. # # # # To name a few: # # - Shamir's Secret sharing scheme # # - Multiplication protocol # # - Integer division # # - Comparison protocol # ##################################################################################################### import numpy as np import matplotlib.pylab as plt import csv import sys #Creates shares of secrets using Shamir's secret sharing scheme. def secretsharing(x, t, n): shares = [] c = [F.random_element() for i in range(t)] for i in range(1, n+1): s = x for j in range(1,t+1): s = F(s) + F(c[j-1]) * F(i)**F(j) shares.append(s) return np.array(shares) #Creates the "recombination"-vector used to reconstruct a secret from its shares. def basispoly(n): r = [] C = range(1,n+1) for i in range(1,n+1): c = [k for k in C if k != i] p = 1 for j in range(n-1): p *= -F(c[j]) / (F(i)-F(c[j])) r.append(F(p)); return r # reconstruct secret. def rec(x): res = 0 n = len(x) y = basispoly(n) for i in range(len(x)): res += F(x[i] * y[i]) return res #Create a Beaver's triplet. def triplet(n, t): r = basispoly(n) matrixA = np.zeros((n,n)) matrixB = np.zeros((n,n)) for i in range(n): a = F.random_element() b = F.random_element() matrixA[:,i] = secretsharing( a, t, n) matrixB[:,i] = secretsharing(b, t, n) sharesA = np.sum(matrixA, 1) sharesA = [F(i) for i in sharesA] sharesB = np.sum(matrixB, 1) sharesB = [F(i) for i in sharesB] sharesC = mult(sharesA,sharesB) #[F(i*j) for (i,j) in zip(sharesA, sharesB)] return [sharesA, sharesB, sharesC] #Secure mulitplication of two secrets. -Resulting polynomial is still degree t. def mult(a,b): r = basispoly(n) h = np.array(a)*np.array(b) hs = [] for i in range(n): hs.append(secretsharing(h[i],t,n)) s = secretsharing(0,t,n) for i in range(n): s+= np.array(r[i]) * np.array(hs[i]) return np.array(s) #Integer division with public constant #def divPub(x, k, l): # Should return a random bit unknown to all parties. Does not seem to work for huge field. def randomBit(n, t): rsq = 0 while rsq == 0: #If r**2 = 0, the parties retry r_shares = triplet(n,t)[0] #The parties compute an unknown random number, r, and computes r**2 print rec(r_shares) rsq_shares = mult(r_shares,r_shares) print rec(rsq_shares) #salpha, sbeta, sgamma = triplet(n, t) #d = dot(F, basis, [F(i-j) for (i,j) in zip(r_shares, salpha)]) #e = dot(F, basis, [F(i-j) for (i,j) in zip(r_shares, sbeta)]) #rsq_shares = [F(d*e + d*i + e*j + k) for (i,j,k) in zip(sbeta, salpha, sgamma)] rsq = rec(rsq_shares) #r**2 is revealed r_prime = sqrt(rsq) print 'rp', r_prime #r = np.random.choice(2,1)[0] return (F(1/2) * (r_prime**(-1) * np.array(r_shares) + 1)) #Return both sharings of a random number and the bitwise sharing of the same random number. def randomBitNumber(F, n, t, l): sharings = [] for i in range(l): sharings.append(randomBit(F,n,t)) sharings_reverse = sharings[::-1] r_shares = [] for i in range(n): s = 0 for j in range(l): s += 2**j * rec(sharings_reverse[j]) r_shares.append(s) return sharings, r_shares #Computing the inverse of a field element. def inverse(x, r): w_s = [i*j for i,j in zip(x,r)] w = rec(w_s) return w**(-1) * np.array(r) #Bitwise AND operation def AND(a): p = a[0] for i in range(1,len(a)): p = mult(p, a[i]) return p #Bitwise XOR operation def xor(a1,b1): d = np.array(a1) - np.array(b1) return mult(d,d) #Bitwise OR operation def OR(a): p = 1 - a[0] for i in range(1,len(a)): p = mult(p, (1 - a[i])) return 1 - p #Converting bitwise sharing's of a secret to ordinary sharings. def dec(a): s = 0 for i in range(len(a)): s+= a[i] * 2**i return s #Computing the truth of whether a < b. def lessthan(a,b): e = [] for i in range(len(a)): e.append(xor(a[i],b[i])) f = [] for i in range(len(e)): f.append(OR(e[0:i+1])) g = [f[0]] for i in range(0,len(f)-1): g.append(f[i+1] - f[i]) h = [] for i in range(len(g)): h.append(mult(g[i], b[i])) y = [] for i in range(len(h[0])): s = 0 for j in range(len(h)): s += h[j][i] y.append(s) return y #Dividing a secret with a public integer def down(F, q, k,t,n): #while True: # r = F.random_element() # if r < 10*dot(F,basis,q): # break r = F(900)# r_T = np.trunc(np.array(r,dtype = float) *k**(-1)) r_s = secretsharing(r,t,n) r_Ts = secretsharing(F(r_T),t,n) z_s = q + r_s z = rec(z_s) z_T = np.trunc(np.array(z, dtype = float) * k**(-1)) z_Ts = secretsharing(F(z_T),t,n) return z_Ts - r_Ts #Computing PRE-OR more efficiently. def pre(F, d, r, rb, n, t, l): ### ADD RANDOMNESS basis = basispoly(n) rb.insert(0, list(secretsharing(0,t,n))) es = np.array(d) + np.array(r) e = rec(es) ebit = list(bin(e)[2:].rjust(l+1,'0')) ebs = [] for b in ebit: ebs.append(secretsharing(b,t,n)) ### Create H ep = [xor(i,j) for i,j in zip(ebs,rb)] E = [] for i in range(len(ep)): E.append(OR(ep[0:i+1])) h = [i - mult(i,j) for i,j in zip (rb, ebs)] hp = [] for i in range(len(h)): if i ==0: Es = 0 else: Es = E[i-1] hp.append(h[i] + 1 - Es) H = [] for i in range(len(hp)): H.append(AND(hp[0:i+1])) ### Calc Hld Hp = [] for i in range(len(H)): Hp.append(xor(H[i],1)) e0 = [] r0 = [] for i in range(len(Hp)): e0.append(mult(Hp[i],ebs[i])) r0.append(mult(Hp[i], rb[i])) B1 = lessthan(e0,r0) H_aim = [1-i for i in B1] ### Trans H into D D_ = [] H2 = H[1:] H2.append(secretsharing(0,t,n)) N = [1-i for i in H_aim] for i in range(len(H)): D_.append(mult(H_aim, H2[i]) + mult(N,H[i])) D = [1-i for i in D_] D = D[::-1] y = [] for i in range(len(D[0])): b = [] for j in range(len(D)): b.append(D[j][i]) y.append(dec(b)) return y #Integer division def integerDiv(F, d, nom, n, t, l, k): r1,r2,c = triplet(n, t) rr = 0 while rr == 0: rbits, r = randomBitNumber(F, n, t, l) rr = rec(r) D = pre(F, d, r, rbits, n, t, l) shift = [1+i for i in D] shiftInv = inverse(shift, r) a1 = [i-j for i,j in zip(shift, d)] p = mult(a1,shiftInv) sum_s = [] prod = secretsharing(1,t,n) s = np.array(n*[0]) for i in range(1,l): prod = mult(prod, p) s = s + np.array(prod) s = [1+i for i in s] a_hat = [2**k*i for i in shiftInv] a_hat = mult(a_hat,s) q_hat = mult(nom,a_hat) q = down(F, q_hat, 2**k,t,n) r = nom - mult(d,q) dd = rec(d+d) rd = rec(r+d) ep = int(dd <= rd) em = int(rec(d) > rd) y = q + ep - em return y ###############################MAIN ########################################## ##### CONSTANTS ################################ #set_random_seed(1) #np.random.seed(1) m = random_prime(2**250, None, 2**247) # Prime number, large enough that the probability of chossing a certain number is # negeligable F = GF(m) # Finite field with m elements n = 4 # Number of parties t = 1 # Degree of polynomial l = 13 # bit length of some secrets r = basispoly(n) # Lagrange basispolynomial coefficients print rec(randomBit(n,t))
 625198650575593907948139785385919338907056866949899424242338243422238322\ 81 634581691528657710873525423649129981338764148962396127626508371755731539\ 231 rp sqrt63458169152865771087352542364912998133876414896239612762650837175573\ 1539231 Traceback (click to the left of this block for traceback) ... ValueError: element is not in the prime field 62519865057559390794813978538591933890705686694989942424233824342223832281 634581691528657710873525423649129981338764148962396127626508371755731539231 rp sqrt634581691528657710873525423649129981338764148962396127626508371755731539231 Traceback (most recent call last): # - Multiplication protocol # File "", line 1, in File "/tmp/tmpxi0y4I/___code___.py", line 316, in exec compile(u'print rec(randomBit(n,t)) File "", line 1, in File "/tmp/tmpxi0y4I/___code___.py", line 46, in rec res += F(x[i] * y[i]) File "sage/structure/parent.pyx", line 1097, in sage.structure.parent.Parent.__call__ (/usr/local/sage-6.9/src/build/cythonized/sage/structure/parent.c:9873) File "sage/rings/finite_rings/hom_prime_finite_field.pyx", line 44, in sage.rings.finite_rings.hom_prime_finite_field.SectionFiniteFieldHomomorphism_prime._call_ (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/hom_prime_finite_field.c:3132) File "/usr/local/sage-6.9/local/lib/python2.7/site-packages/sage/rings/finite_rings/integer_mod_ring.py", line 1167, in _element_constructor_ return integer_mod.IntegerMod(self, x) File "sage/rings/finite_rings/integer_mod.pyx", line 182, in sage.rings.finite_rings.integer_mod.IntegerMod (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/integer_mod.c:4160) File "sage/rings/finite_rings/integer_mod.pyx", line 1716, in sage.rings.finite_rings.integer_mod.IntegerMod_gmp.__init__ (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/integer_mod.c:20448) File "sage/structure/parent.pyx", line 1097, in sage.structure.parent.Parent.__call__ (/usr/local/sage-6.9/src/build/cythonized/sage/structure/parent.c:9873) File "sage/structure/coerce_maps.pyx", line 237, in sage.structure.coerce_maps.NamedConvertMap._call_ (/usr/local/sage-6.9/src/build/cythonized/sage/structure/coerce_maps.c:5944) File "sage/rings/finite_rings/element_pari_ffelt.pyx", line 937, in sage.rings.finite_rings.element_pari_ffelt.FiniteFieldElement_pari_ffelt._integer_ (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/element_pari_ffelt.c:9355) File "sage/rings/finite_rings/element_pari_ffelt.pyx", line 920, in sage.rings.finite_rings.element_pari_ffelt.FiniteFieldElement_pari_ffelt.lift (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/element_pari_ffelt.c:9228) ValueError: element is not in the prime field
random_prime(2**250, None, 2**247)
 696205902887547871652159517929547663282608944037200206897421778089495766\ 033 696205902887547871652159517929547663282608944037200206897421778089495766033
import numpy as np import matplotlib.pylab as plt import csv import sys #Creates shares of secrets using Shamir's secret sharing scheme. def secretsharing(x, t, n): shares = [] c = [F.random_element() for i in range(t)] for i in range(1, n+1): s = x for j in range(1,t+1): s = F(s) + F(c[j-1]) * F(i)**F(j) shares.append(s) return np.array(shares) g #Creates the "recombination"-vector used to reconstruct a secret from its shares. def basispoly(n): r = [] C = range(1,n+1) for i in range(1,n+1): c = [k for k in C if k != i] p = 1 for j in range(n-1): p *= -F(c[j]) / (F(i)-F(c[j])) r.append(F(p)); return r # reconstruct secret. def rec(x): res = 0 n = len(x) y = basispoly(n) for i in range(len(x)): res += F(x[i] * y[i]) return res #Create a Beaver's triplet. def triplet(n, t): r = basispoly(n) matrixA = np.zeros((n,n)) matrixB = np.zeros((n,n)) for i in range(n): a = F.random_element() b = F.random_element() matrixA[:,i] = secretsharing( a, t, n) matrixB[:,i] = secretsharing(b, t, n) sharesA = np.sum(matrixA, 1) sharesA = [F(i) for i in sharesA] sharesB = np.sum(matrixB, 1) sharesB = [F(i) for i in sharesB] sharesC = mult(sharesA,sharesB) #[F(i*j) for (i,j) in zip(sharesA, sharesB)] return [sharesA, sharesB, sharesC] #Secure mulitplication of two secrets. -Resulting polynomial is still degree t. def mult(a,b): r = basispoly(n) h = np.array(a)*np.array(b) hs = [] for i in range(n): hs.append(secretsharing(h[i],t,n)) s = secretsharing(0,t,n) for i in range(n): s+= np.array(r[i]) * np.array(hs[i]) return np.array(s) #Integer division with public constant #def divPub(x, k, l): m = random_prime(2**250, None, 2**247) # Prime number, large enough that the probability of chossing a certain number is # negeligable F = GF(m) # Finite field with m elements n = 4 # Number of parties t = 1 # Degree of polynomial l = 13 # bit length of some secrets r = basispoly(n) # Lagrange basispolynomial coefficients def randomBit(n, t): rsq = 0 while rsq == 0: #If r**2 = 0, the parties retry r_shares = triplet(n,t)[0] #The parties compute an unknown random number, r, and computes r**2 rsq_shares = mult(r_shares,r_shares) rsq = rec(rsq_shares) #r**2 is revealed r_prime = sqrt(rsq) return (F(1/2) * (r_prime**(-1) * np.array(r_shares) + 1)) print rec(randomBit(n,t))
 Traceback (click to the left of this block for traceback) ... ValueError: element is not in the prime field Traceback (most recent call last): shares = [] File "", line 1, in File "/tmp/tmpbir8B8/___code___.py", line 97, in exec compile(u'print rec(randomBit(n,t)) File "", line 1, in File "/tmp/tmpbir8B8/___code___.py", line 36, in rec res += F(x[i] * y[i]) File "sage/structure/parent.pyx", line 1097, in sage.structure.parent.Parent.__call__ (/usr/local/sage-6.9/src/build/cythonized/sage/structure/parent.c:9873) File "sage/rings/finite_rings/hom_prime_finite_field.pyx", line 44, in sage.rings.finite_rings.hom_prime_finite_field.SectionFiniteFieldHomomorphism_prime._call_ (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/hom_prime_finite_field.c:3132) File "/usr/local/sage-6.9/local/lib/python2.7/site-packages/sage/rings/finite_rings/integer_mod_ring.py", line 1167, in _element_constructor_ return integer_mod.IntegerMod(self, x) File "sage/rings/finite_rings/integer_mod.pyx", line 182, in sage.rings.finite_rings.integer_mod.IntegerMod (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/integer_mod.c:4160) File "sage/rings/finite_rings/integer_mod.pyx", line 1716, in sage.rings.finite_rings.integer_mod.IntegerMod_gmp.__init__ (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/integer_mod.c:20448) File "sage/structure/parent.pyx", line 1097, in sage.structure.parent.Parent.__call__ (/usr/local/sage-6.9/src/build/cythonized/sage/structure/parent.c:9873) File "sage/structure/coerce_maps.pyx", line 237, in sage.structure.coerce_maps.NamedConvertMap._call_ (/usr/local/sage-6.9/src/build/cythonized/sage/structure/coerce_maps.c:5944) File "sage/rings/finite_rings/element_pari_ffelt.pyx", line 937, in sage.rings.finite_rings.element_pari_ffelt.FiniteFieldElement_pari_ffelt._integer_ (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/element_pari_ffelt.c:9355) File "sage/rings/finite_rings/element_pari_ffelt.pyx", line 920, in sage.rings.finite_rings.element_pari_ffelt.FiniteFieldElement_pari_ffelt.lift (/usr/local/sage-6.9/src/build/cythonized/sage/rings/finite_rings/element_pari_ffelt.c:9228) ValueError: element is not in the prime field

b########################################################################### #This cells defines functions that are more specifically used to # #compute the recursive least squares equations securely. # ########################################################################### #Should and could, be implemented securely. def less(a,b): return int( rec(a) < rec(b) ) def reVector(b, n): Nvector = [] for i in range(2): liste = [] for j in range(n): liste.append(b[j][i][0]) Nvector.append(F(rec(liste)))#[b[0][1][i][0],b[1][1][i][0],b[2][1][i][0],b[3][1][i][0]])) return matrix(2,1,Nvector) def ILS(F, P, x, w, y, C, t, n, i): g = matVec(P,x) e0 = vecVec(x,w) e0 = down(F,e0,C,t,n) e = [k-j for k,j in zip(y, e0)] a = [mult(k,e) for k in g] w = [i+j for i,j in zip(w,a)] return w def estP(F,P, x, C,t,n): xP = vecMat(x,P) a1 = vecVec(xP, x) denom = C + a1 Px = matVec(P,x) Part = VecVec(Px, xP) return denom, Part #Computing securely with coorect degree polynomial: x.tranpose * y def vecVec(x,y): p = len(x) s = secretsharing(0,t,n) for i in range(p): s+= mult(x[i], y[i]) return s #Computing securely with correct degree polynomial: x.transpose * A def vecMat(x,A): p = len(x) res = [] for i in range(p): s = secretsharing(0,t,n) for j in range(p): s+= mult(x[j], A[i+j*p]) res.append(s) return res #Computing securely with correct degree polynomial: A * x def matVec(A,x): res = [] p = len(x) for i in range(p): s = secretsharing(0,t,n) for j in range(p): s += mult(x[j], A[j+(i*p)]) res.append(s) return res #Computing securely with correct degree polynomial: x * x.transpose def VecVec(x,y): res = [] for i in range(len(y)): for j in range(len(x)): res.append(mult(y[i],x[j])) return res
########################################################################################### # This cells computes SECURELY the recursice least square equations using the functions # # from cell 1 and 2 (thus cell 1 and 2 must run before attempting to run this cell). # ########################################################################################### np.random.seed(2) obs = 70 n = 4 # Numbers of parties t = 1 # Corrupted parties p = 2 # Dimension (for one dimensional system, p = 2) C = 2**7 # Scaling constant P = C*np.eye(p) # Initialize P Ps = [] # List for holding shares of P l = 13 # Maximum number of bits for any secret nominator and denominator k=l**2+l # Scaling for 1/d #listObs = [obs*[1]] # 1 dimensional system listObs = [] # Multi dimensional system #params = [2] # 1 dimensional system params = [] # Multi dimensional system for j in range(p): # For 1 dimensional system write p-1. listObs.append(np.random.randint(0,5,obs)) #Observations of input params.append(np.random.randint(1,5,1)) xobs = np.array(listObs).transpose() params = np.array(params).flatten() print 'True parameters: ', params print pretty_print(html('<u> Time | Estimate </u>')) y = [sum(i * params) for i in xobs] #Observations of output noice = np.random.normal(0,2,obs) #Noice y = y+noice #Add noice y = np.array(y,dtype = int) #Truncate Y #The P matrix, will be a list of the secret entries. --> Ps = [P[1,1], P[1,2], P[2,1], P[2,2]]. for i in range(p): for j in range(p): Ps.append(secretsharing(P[i][j],t,n)) w = [secretsharing(0,t,n) for i in range(p)] #Initialization of the parameter estimate wOPEN = [] for i in range(obs): x = xobs[i] #Observation vector for 3 X = [list(secretsharing(j,t,n)) for j in x] #Secret sharing of observations denom, Ppart = estP(F,Ps, X,C,t,n) # TJECK FOR WRAP-AROUND ZERO (COULD and SHOULD be implemented securely) b = [less(Ppart[j]*-1,Ppart[j]) for j in range(p**2)] #REMOVE SIGN for j in range(p**2): if b[j]: Ppart[j] = -Ppart[j] #PERFORM DIVISION P2 = [integerDiv(F, denom, np.array(j), n, t, l, k) for j in Ppart] #ADD SIGN for j in range(p**2): if b[j]: P2[j] = -P2[j] Ps = [h-j for h,j in zip(Ps, P2)] y2 = secretsharing(y[i],t,n) w = ILS(F, Ps, X, w, y2, C,t,n,i) w2 = np.copy(w) w_open = [] #TJECK WRAP-AROUND ZERO b2 = [less(j*-1, j) for j in w2] for j in range(p): if b2[j]: w2[j] = -w2[j] #REVEAL (FLOATING POINT) RESULT WITH ADDED SIGN for j in w2: w_open.append(np.array(rec(j),dtype = float)/C) for j in range(p): if b2[j]: w_open[j] = -w_open[j] w_open = np.array(w_open) print ' {:>2} | {} '.format(i, w_open) #print i'th estimate. wOPEN.append(w_open) print '---------------------------------------' e = [] for j in range(obs): e.append(sum((params - wOPEN[j])**2) * 1./p) #print '({},{:.10f})'.format(j, e[j]) #PLOTTING PURPOSE list_plot(e)
 True parameters: [3 2] Time | Estimate 0 | [ 0. 2.] 1 | [ 0. 2.] 2 | [ 2.7421875 2.2109375] 3 | [ 2.8984375 2.3671875] 4 | [ 2.5234375 2.4609375] 5 | [ 2.5234375 2.4609375] 6 | [ 2.5234375 2.4609375] 7 | [ 2.5703125 2.3671875] 8 | [ 2.921875 2.3671875] 9 | [ 2.921875 2.3671875] 10 | [ 3.390625 2.0859375] 11 | [ 3.3046875 2.0625 ] 12 | [ 3.21875 2.0390625] 13 | [ 3.40625 1.8984375] 14 | [ 3.21875 2.0234375] 15 | [ 3.234375 2.0703125] 16 | [ 3.234375 2.0703125] 17 | [ 3.0703125 2.09375 ] 18 | [ 3.0390625 2. ] 19 | [ 3.0625 1.921875] 20 | [ 3.03125 2.109375] 21 | [ 2.984375 2.078125] 22 | [ 2.9609375 2.03125 ] 23 | [ 2.9609375 2.03125 ] 24 | [ 2.8671875 1.84375 ] 25 | [ 2.9453125 1.875 ] 26 | [ 2.921875 1.953125] 27 | [ 3.0625 1.859375] 28 | [ 3.0625 1.859375] 29 | [ 3.125 1.859375] 30 | [ 3.125 1.859375] 31 | [ 3.140625 1.765625] 32 | [ 3.140625 1.765625] 33 | [ 3.140625 1.765625] 34 | [ 3.1484375 1.78125 ] 35 | [ 3.1328125 1.8125 ] 36 | [ 3.0078125 1.8125 ] 37 | [ 2.9921875 1.78125 ] 38 | [ 2.9921875 1.78125 ] 39 | [ 2.8515625 2.0625 ] 40 | [ 2.8984375 2.03125 ] 41 | [ 2.8828125 2.0625 ] 42 | [ 3.0234375 1.96875 ] 43 | [ 3.0625 1.890625] 44 | [ 3.0625 1.890625] 45 | [ 3.0625 1.890625] 46 | [ 3.15625 1.90625] 47 | [ 3.015625 2. ] 48 | [ 3.0859375 1.953125 ] 49 | [ 3.0390625 2.0234375] 50 | [ 3.1171875 2.0234375] 51 | [ 3.1796875 1.9296875] 52 | [ 3.1328125 1.8828125] 53 | [ 3.1015625 1.890625 ] 54 | [ 3.0546875 1.9609375] 55 | [ 3.1171875 2.0234375] 56 | [ 3.1171875 2.1015625] 57 | [ 3.1328125 2.078125 ] 58 | [ 3.2265625 2.015625 ] 59 | [ 3.3203125 2.03125 ] 60 | [ 3.296875 2.046875] 61 | [ 3.34375 2.09375] 62 | [ 3.359375 2.0703125] 63 | [ 3.359375 2.1484375] 64 | [ 3.453125 1.9296875] 65 | [ 3.171875 2.1171875] 66 | [ 3.234375 2.0234375] 67 | [ 3.15625 2.1796875] 68 | [ 3.1484375 2.2109375] 69 | [ 3.1640625 2.2265625] --------------------------------------- True parameters: [3 2] Time | Estimate 0 | [ 0. 2.] 1 | [ 0. 2.] 2 | [ 2.7421875 2.2109375] 3 | [ 2.8984375 2.3671875] 4 | [ 2.5234375 2.4609375] 5 | [ 2.5234375 2.4609375] 6 | [ 2.5234375 2.4609375] 7 | [ 2.5703125 2.3671875] 8 | [ 2.921875 2.3671875] 9 | [ 2.921875 2.3671875] 10 | [ 3.390625 2.0859375] 11 | [ 3.3046875 2.0625 ] 12 | [ 3.21875 2.0390625] 13 | [ 3.40625 1.8984375] 14 | [ 3.21875 2.0234375] 15 | [ 3.234375 2.0703125] 16 | [ 3.234375 2.0703125] 17 | [ 3.0703125 2.09375 ] 18 | [ 3.0390625 2. ] 19 | [ 3.0625 1.921875] 20 | [ 3.03125 2.109375] 21 | [ 2.984375 2.078125] 22 | [ 2.9609375 2.03125 ] 23 | [ 2.9609375 2.03125 ] 24 | [ 2.8671875 1.84375 ] 25 | [ 2.9453125 1.875 ] 26 | [ 2.921875 1.953125] 27 | [ 3.0625 1.859375] 28 | [ 3.0625 1.859375] 29 | [ 3.125 1.859375] 30 | [ 3.125 1.859375] 31 | [ 3.140625 1.765625] 32 | [ 3.140625 1.765625] 33 | [ 3.140625 1.765625] 34 | [ 3.1484375 1.78125 ] 35 | [ 3.1328125 1.8125 ] 36 | [ 3.0078125 1.8125 ] 37 | [ 2.9921875 1.78125 ] 38 | [ 2.9921875 1.78125 ] 39 | [ 2.8515625 2.0625 ] 40 | [ 2.8984375 2.03125 ] 41 | [ 2.8828125 2.0625 ] 42 | [ 3.0234375 1.96875 ] 43 | [ 3.0625 1.890625] 44 | [ 3.0625 1.890625] 45 | [ 3.0625 1.890625] 46 | [ 3.15625 1.90625] 47 | [ 3.015625 2. ] 48 | [ 3.0859375 1.953125 ] 49 | [ 3.0390625 2.0234375] 50 | [ 3.1171875 2.0234375] 51 | [ 3.1796875 1.9296875] 52 | [ 3.1328125 1.8828125] 53 | [ 3.1015625 1.890625 ] 54 | [ 3.0546875 1.9609375] 55 | [ 3.1171875 2.0234375] 56 | [ 3.1171875 2.1015625] 57 | [ 3.1328125 2.078125 ] 58 | [ 3.2265625 2.015625 ] 59 | [ 3.3203125 2.03125 ] 60 | [ 3.296875 2.046875] 61 | [ 3.34375 2.09375] 62 | [ 3.359375 2.0703125] 63 | [ 3.359375 2.1484375] 64 | [ 3.453125 1.9296875] 65 | [ 3.171875 2.1171875] 66 | [ 3.234375 2.0234375] 67 | [ 3.15625 2.1796875] 68 | [ 3.1484375 2.2109375] 69 | [ 3.1640625 2.2265625] ---------------------------------------
########################################################################################### # This cells computes NON-SECURELY the recursice least square equations. # # Thus FLOATING-point numbers are used in these equations. Can be used as a comparison # # for the results from cell 3. # ########################################################################################### import numpy as np np.random.seed(2) obs = 70 p = 2 #listObs = [obs*[1]] #1 dimensional system listObs = [] #Multi dimensional system #params = [2] #1 dimensional system params = [] #Multi dimensional system for j in range(p): listObs.append(np.random.randint(0,5,obs)) #Observations of input params.append(np.random.randint(1,5,1)) xobs = np.array(listObs).transpose() params = np.array(params).flatten() print 'True parameters: ', params print pretty_print(html('<u> Time | Estimate </u>')) y = [sum(i * params) for i in xobs] #Observations of output noice = np.random.normal(0,2,obs) #Noice y = y+noice #Add noice P = matrix.identity(p) w = matrix(p,1,p*[0]) wOPEN = [] for i in range(obs): xdut = [int(j) for j in xobs[i]] #x = matrix(p,1,[1, x1[i]]) #Observation vector for 1 dimension x = matrix(p,1, xdut) #Observation vector for multi dimensions denom = (1 + x.transpose()* P * x) P = P - matrix(denom[0][0]**(-1) *( P * x * x.transpose() * P)) g = P * x e = y[i] - x.transpose() * w w = matrix(w + g * e) w3 = np.array(w) print ' {:>2} | {} '.format(i, w_open) #print i'th estimate. wOPEN.append(w3) print '---------------------------------------' params = params.reshape(p,1) e = [] for j in range(obs): e.append(sum((params - wOPEN[j])**2) * 1./p) #print '({},{:.10f})'.format(j, e[j]) list_plot(e)
 True parameters: [3 2] Time | Estimate 0 | [ 1.875 2.3671875] 1 | [ 1.875 2.3671875] 2 | [ 1.875 2.3671875] 3 | [ 1.875 2.3671875] 4 | [ 1.875 2.3671875] 5 | [ 1.875 2.3671875] 6 | [ 1.875 2.3671875] 7 | [ 1.875 2.3671875] 8 | [ 1.875 2.3671875] 9 | [ 1.875 2.3671875] 10 | [ 1.875 2.3671875] 11 | [ 1.875 2.3671875] 12 | [ 1.875 2.3671875] 13 | [ 1.875 2.3671875] 14 | [ 1.875 2.3671875] 15 | [ 1.875 2.3671875] 16 | [ 1.875 2.3671875] 17 | [ 1.875 2.3671875] 18 | [ 1.875 2.3671875] 19 | [ 1.875 2.3671875] 20 | [ 1.875 2.3671875] 21 | [ 1.875 2.3671875] 22 | [ 1.875 2.3671875] 23 | [ 1.875 2.3671875] 24 | [ 1.875 2.3671875] 25 | [ 1.875 2.3671875] 26 | [ 1.875 2.3671875] 27 | [ 1.875 2.3671875] 28 | [ 1.875 2.3671875] 29 | [ 1.875 2.3671875] 30 | [ 1.875 2.3671875] 31 | [ 1.875 2.3671875] 32 | [ 1.875 2.3671875] 33 | [ 1.875 2.3671875] 34 | [ 1.875 2.3671875] 35 | [ 1.875 2.3671875] 36 | [ 1.875 2.3671875] 37 | [ 1.875 2.3671875] 38 | [ 1.875 2.3671875] 39 | [ 1.875 2.3671875] 40 | [ 1.875 2.3671875] 41 | [ 1.875 2.3671875] 42 | [ 1.875 2.3671875] 43 | [ 1.875 2.3671875] 44 | [ 1.875 2.3671875] 45 | [ 1.875 2.3671875] 46 | [ 1.875 2.3671875] 47 | [ 1.875 2.3671875] 48 | [ 1.875 2.3671875] 49 | [ 1.875 2.3671875] 50 | [ 1.875 2.3671875] 51 | [ 1.875 2.3671875] 52 | [ 1.875 2.3671875] 53 | [ 1.875 2.3671875] 54 | [ 1.875 2.3671875] 55 | [ 1.875 2.3671875] 56 | [ 1.875 2.3671875] 57 | [ 1.875 2.3671875] 58 | [ 1.875 2.3671875] 59 | [ 1.875 2.3671875] 60 | [ 1.875 2.3671875] 61 | [ 1.875 2.3671875] 62 | [ 1.875 2.3671875] 63 | [ 1.875 2.3671875] 64 | [ 1.875 2.3671875] 65 | [ 1.875 2.3671875] 66 | [ 1.875 2.3671875] 67 | [ 1.875 2.3671875] 68 | [ 1.875 2.3671875] 69 | [ 1.875 2.3671875] --------------------------------------- True parameters: [3 2] Time | Estimate 0 | [ 1.875 2.3671875] 1 | [ 1.875 2.3671875] 2 | [ 1.875 2.3671875] 3 | [ 1.875 2.3671875] 4 | [ 1.875 2.3671875] 5 | [ 1.875 2.3671875] 6 | [ 1.875 2.3671875] 7 | [ 1.875 2.3671875] 8 | [ 1.875 2.3671875] 9 | [ 1.875 2.3671875] 10 | [ 1.875 2.3671875] 11 | [ 1.875 2.3671875] 12 | [ 1.875 2.3671875] 13 | [ 1.875 2.3671875] 14 | [ 1.875 2.3671875] 15 | [ 1.875 2.3671875] 16 | [ 1.875 2.3671875] 17 | [ 1.875 2.3671875] 18 | [ 1.875 2.3671875] 19 | [ 1.875 2.3671875] 20 | [ 1.875 2.3671875] 21 | [ 1.875 2.3671875] 22 | [ 1.875 2.3671875] 23 | [ 1.875 2.3671875] 24 | [ 1.875 2.3671875] 25 | [ 1.875 2.3671875] 26 | [ 1.875 2.3671875] 27 | [ 1.875 2.3671875] 28 | [ 1.875 2.3671875] 29 | [ 1.875 2.3671875] 30 | [ 1.875 2.3671875] 31 | [ 1.875 2.3671875] 32 | [ 1.875 2.3671875] 33 | [ 1.875 2.3671875] 34 | [ 1.875 2.3671875] 35 | [ 1.875 2.3671875] 36 | [ 1.875 2.3671875] 37 | [ 1.875 2.3671875] 38 | [ 1.875 2.3671875] 39 | [ 1.875 2.3671875] 40 | [ 1.875 2.3671875] 41 | [ 1.875 2.3671875] 42 | [ 1.875 2.3671875] 43 | [ 1.875 2.3671875] 44 | [ 1.875 2.3671875] 45 | [ 1.875 2.3671875] 46 | [ 1.875 2.3671875] 47 | [ 1.875 2.3671875] 48 | [ 1.875 2.3671875] 49 | [ 1.875 2.3671875] 50 | [ 1.875 2.3671875] 51 | [ 1.875 2.3671875] 52 | [ 1.875 2.3671875] 53 | [ 1.875 2.3671875] 54 | [ 1.875 2.3671875] 55 | [ 1.875 2.3671875] 56 | [ 1.875 2.3671875] 57 | [ 1.875 2.3671875] 58 | [ 1.875 2.3671875] 59 | [ 1.875 2.3671875] 60 | [ 1.875 2.3671875] 61 | [ 1.875 2.3671875] 62 | [ 1.875 2.3671875] 63 | [ 1.875 2.3671875] 64 | [ 1.875 2.3671875] 65 | [ 1.875 2.3671875] 66 | [ 1.875 2.3671875] 67 | [ 1.875 2.3671875] 68 | [ 1.875 2.3671875] 69 | [ 1.875 2.3671875] ---------------------------------------