Privacy Preserving Homomorphic Linear Regression
Problem Setup
This example, based on the homomorpheR package, was communicated to me by Kee Siong Ng.
FIUholds the response vectoryRE1holds the datax[, 1:3]RE2RE2 holds the datax[,4]
Together, they want to learn a linear model to predict y using x but in such a way that
- nobody sees each other’s data,
RE1holds the coefficients for the variables it holds, which is not visible to othersRE2holds the coefficients for the variables it holds, which is not visible to others- Prediction of a new instance requires the collaboration of
RE1andRE2
FIU first generates a public-private key pair using the Paillier scheme. The public key is used by RE1 and RE2 to encrypt their data and do the computations. The private key is visible only to FIU and is used by the FIU to decrypt data sent from the REs.
A Synthetic Data Set
n <- 50
x1 = 11:(10+n)
x2 = runif(n,5,95)
x3 = rbinom(n,1,0.5)
x = data.frame(x1,x2,x3)
x = scale(x)
x = data.frame(1,x) # the intercept is handled using a column of 1's
sigma = 1.4
eps = rnorm(x1, 0, sigma) # generate noise vector
b = c(17, -2.5, 0.5,-5.2) # the real model coefficient
y = b[1] + b[2] * x[,2] + b[3] * x[,3] + b[4] * x[,4] + scale(eps) # target variableBatch Gradient Descent
This example uses a batch gradient descent algorithm for fitting the model.
lm_sgd <- function(iter, rate) {
theta = c(0,0,0,0)
alpha = rate
for (iter in 1:iter) {
adj = c(0,0,0,0)
for (i in 1:4) {
for (j in 1:n) {
adj[i] = adj[i] + (sum(x[j,] * theta) - y[j]) * x[j,i]
}
# theta[i] = theta[i] - (alpha / n) * adj[i] # adjust as we go is faster
}
for (i in 1:4) theta[i] = theta[i] - (alpha / n) * adj[i]
print(adj)
}
print(theta)
}As a check, we can run 100 iterations to see if it works.
lm_sgd(100, 0.1)## [1] -850.00000 118.25073 -48.15952 255.53613
## [1] -765.00000 106.44608 -40.56685 230.69506
## [1] -688.50000 95.87435 -34.00341 208.27524
## [1] -619.65000 86.40182 -28.33898 188.04018
## [1] -557.68500 77.90971 -23.45922 169.77652
## [1] -501.91650 70.29246 -19.26379 153.29177
## [1] -451.72485 63.45623 -15.66464 138.41227
## [1] -406.55236 57.31756 -12.58456 124.98134
## [1] -365.897128 51.802222 -9.955889 112.857635
## [1] -329.307416 46.844136 -7.719341 101.913607
## [1] -296.376674 42.384486 -5.823009 92.034188
## [1] -266.739007 38.370887 -4.221473 83.115553
## [1] -240.065106 34.756662 -2.875014 75.064024
## [1] -216.058595 31.500205 -1.748921 67.795075
## [1] -194.4527359 28.5644164 -0.8128824 61.2324350
## [1] -175.00746228 25.91619578 -0.04044792 55.30728536
## [1] -157.5067161 23.5260007 0.5914422 49.9575257
## [1] -141.756044 21.367451 1.102875 45.127119
## [1] -127.580440 19.416979 1.511329 40.765499
## [1] -114.82240 17.65352 1.83200 36.82703
## [1] -103.340156 16.058238 2.078083 33.270540
## [1] -93.006141 14.614279 2.261018 30.058861
## [1] -83.705527 13.306556 2.390715 27.158457
## [1] -75.334974 12.121558 2.475744 24.539060
## [1] -67.801477 11.047176 2.523503 22.173354
## [1] -61.021329 10.072556 2.540369 20.036686
## [1] -54.919196 9.187959 2.531827 18.106804
## [1] -49.427276 8.384644 2.502583 16.363624
## [1] -44.484549 7.654762 2.456670 14.789019
## [1] -40.03609 6.99126 2.39753 13.36663
## [1] -36.032485 6.387795 2.328094 12.081678
## [1] -32.429236 5.838664 2.250848 10.920840
## [1] -29.186312 5.338731 2.167895 9.872079
## [1] -26.267681 4.883371 2.081003 8.924531
## [1] -23.640913 4.468418 1.991654 8.068388
## [1] -21.276822 4.090112 1.901081 7.294797
## [1] -19.149140 3.745065 1.810302 6.595763
## [1] -17.234226 3.430215 1.720152 5.964070
## [1] -15.510803 3.142798 1.631309 5.393203
## [1] -13.959723 2.880316 1.544314 4.877278
## [1] -12.563751 2.640510 1.459595 4.410983
## [1] -11.307375 2.421334 1.377482 3.989522
## [1] -10.176638 2.220940 1.298221 3.608564
## [1] -9.158974 2.037649 1.221988 3.264197
## [1] -8.243077 1.869943 1.148902 2.952890
## [1] -7.418769 1.716444 1.079031 2.671454
## [1] -6.676892 1.575901 1.012403 2.417006
## [1] -6.0092029 1.4471792 0.9490125 2.1869464
## [1] -5.4082826 1.3292482 0.8888269 1.9789246
## [1] -4.8674544 1.2211708 0.8317909 1.7908184
## [1] -4.3807089 1.1220949 0.7778318 1.6207107
## [1] -3.9426380 1.0312455 0.7268632 1.4668697
## [1] -3.5483742 0.9479173 0.6787881 1.3277312
## [1] -3.1935368 0.8714678 0.6335017 1.2018819
## [1] -2.8741831 0.8013118 0.5908941 1.0880452
## [1] -2.5867648 0.7369158 0.5508516 0.9850675
## [1] -2.3280883 0.6777935 0.5132589 0.8919069
## [1] -2.0952795 0.6235009 0.4780000 0.8076216
## [1] -1.8857515 0.5736331 0.4449597 0.7313608
## [1] -1.6971764 0.5278201 0.4140243 0.6623556
## [1] -1.5274588 0.4857242 0.3850822 0.5999113
## [1] -1.3747129 0.4470366 0.3580247 0.5434000
## [1] -1.2372416 0.4114750 0.3327464 0.4922543
## [1] -1.1135174 0.3787814 0.3091455 0.4459611
## [1] -1.0021657 0.3487196 0.2871238 0.4040570
## [1] -0.9019491 0.3210735 0.2665872 0.3661229
## [1] -0.8117542 0.2956452 0.2474459 0.3317800
## [1] -0.7305788 0.2722535 0.2296140 0.3006858
## [1] -0.6575209 0.2507324 0.2130098 0.2725309
## [1] -0.5917688 0.2309298 0.1975557 0.2470353
## [1] -0.5325919 0.2127064 0.1831780 0.2239459
## [1] -0.4793327 0.1959341 0.1698071 0.2030339
## [1] -0.4313995 0.1804958 0.1573774 0.1840924
## [1] -0.3882595 0.1662839 0.1458265 0.1669341
## [1] -0.3494336 0.1531997 0.1350962 0.1513900
## [1] -0.3144902 0.1411526 0.1251312 0.1373069
## [1] -0.2830412 0.1300593 0.1158799 0.1245464
## [1] -0.2547371 0.1198436 0.1072936 0.1129831
## [1] -0.22926336 0.11043512 0.09932677 0.10250396
## [1] -0.20633703 0.10176956 0.09193665 0.09300629
## [1] -0.18570333 0.09378764 0.08508322 0.08439743
## [1] -0.16713299 0.08643494 0.07872903 0.07659346
## [1] -0.15041969 0.07966142 0.07283905 0.06951846
## [1] -0.13537772 0.07342108 0.06738056 0.06310371
## [1] -0.12183995 0.06767163 0.06232300 0.05728705
## [1] -0.10965596 0.06237417 0.05763786 0.05201218
## [1] -0.09869036 0.05749292 0.05329854 0.04722818
## [1] -0.08882132 0.05299496 0.04928024 0.04288892
## [1] -0.07993919 0.04885001 0.04555986 0.03895267
## [1] -0.07194527 0.04503021 0.04211588 0.03538164
## [1] -0.06475075 0.04150991 0.03892827 0.03214160
## [1] -0.05827567 0.03826550 0.03597840 0.02920157
## [1] -0.05244810 0.03527528 0.03324894 0.02653350
## [1] -0.04720329 0.03251922 0.03072377 0.02411196
## [1] -0.04248296 0.02997891 0.02838790 0.02191395
## [1] -0.03823467 0.02763741 0.02622744 0.01991862
## [1] -0.03441120 0.02547910 0.02422945 0.01810708
## [1] -0.03097008 0.02348958 0.02238193 0.01646223
## [1] -0.02787307 0.02165563 0.02067374 0.01496856
## [1] -0.02508577 0.01996504 0.01909455 0.01361203
## [1] 16.9995485 -2.4880167 0.3689761 -5.2744076Some Tools
Some tools are needed for performing computations using real numbers. Encryption schemes all work on large integers.
library(homomorpheR)
keypair = PaillierKeyPair$new(modulusBits = 1024)
pubkey = keypair$pubkey
privkey = keypair$getPrivateKey()
# Convenience function
zipf = function(x1,x2,x3) {
res = c()
for (i in 1:length(x1)) res = c(res,x1[i],x2[i],x3[i])
res
}
# zipf(5,10,1)
# zipf(c(1,2,3), c(4,5,6), c(1,1,1))
######################################################################
# Paillier is defined only on positive integers.
# We approximate a floating point number X by 3 integers (a,b,c)
# such that X is approximately equal to (a + b/largenum) / 10^c
#
# Example: Assuming largenum = 2^10, 5.2 can be represented by
# (5, 204, 0) = (5 + 204/1024) / 1 = 5.19922
# or (50, 2040, 1) = (50 + 2040/1024) / 10 = 5.19922
#
# In general, the larger largenum is, the better the approximation.
#################
largenum = 2^50
encodeFloat = function(x, e=0) {
x = x * 10^e
x1 = floor(x)
x2 = x - x1
x2 = floor(x2 * largenum)
exps = rep(e, length(x))
zipf(x1,x2,exps)
}
# encodeFloat(5.2)
# encodeFloat(5.2, 2)
# encodeFloat(c(4.5,5.2,7), 1)
# Some useful constants
zero = pubkey$encrypt('0')
one = pubkey$encrypt('1')
############################################################################
# We encrypt each number (a,b,c) by encrypting a and b with the public key
# The input can be a vector of numbers
################
encrypt = function(z, e=0) {
y = encodeFloat(z, e)
res = rep(zero, length(y))
for (i in 1:length(y)) {
if (i %% 3 == 0) { res[i] = y[i] }
else { res[i] = pubkey$encrypt(as.character(y[i])) }
}
res
}
############################################################################
# We decrypt each encrypted number (enc(a), enc(b),c) using the private key
# by computing (dec(enc(a)) + dec(enc(b))/largenum) / 10^c
# The input can be a vector of numbers.
# We need to deal with negative numbers too.
##################
decrypt = function(en) {
uen = en
for (i in 1:length(uen)) {
if (i %% 3 != 0) uen[i] = privkey$decrypt(en[i])
}
# deal with negative results;
# see https://crypto.stackexchange.com/questions/19457/how-can-i-do-minus-on-plaintexts-in-the-paillier-cryptosystem
for (i in 1:length(uen)) {
if (i %% 3 != 0 && # don't process exponents
uen[i] >= as.double(pubkey$n/3.0)) { uen[i] = uen[i] - pubkey$n }
}
res = c()
for (i in seq(1, length(uen), 3)) {
res = c(res,
(as.double(uen[i]) + as.double(uen[i+1]/(largenum))) / as.double(10^uen[i+2]))
}
return(res)
}
# decrypt(encrypt(4.7))
# decrypt(encrypt(0.5))
# decrypt(encrypt(-2.3))
# (z = encrypt(c(5.2,4.7,-2.6)))
# decrypt(z)
############################################################################
# We next define addition of two encrypted numbers: (enc(a1), enc(b1), c1) and (enc(a2), enc(b2), c2).
# Addition can be done on the encrypted a's and encrypted b's if c1 = c2. We make c1 = c2 true by
# multiplying by powers of 10 when necessary.
# The need to make c1 = c2 is the reason why we can't have c1 and c2 encrypted.
# The function works on vectors.
###############
addenc = function(x, y) {
res = x
for (i in seq(1, length(x), 3)) {
if (x[i+2] != y[i+2]) { # if different exponent, make the exponents equal by multiplying
res[i+2] = max(x[i+2], y[i+2])
xdiff = res[i+2] - x[i+2]
ydiff = res[i+2] - y[i+2]
if (xdiff > 0) {
x[i] = pubkey$mult(x[i], 10^(xdiff))
x[i+1] = pubkey$mult(x[i+1], 10^(xdiff))
}
if (ydiff > 0) {
y[i] = pubkey$mult(y[i], 10^(ydiff))
y[i+1] = pubkey$mult(y[i+1], 10^(ydiff))
}
}
res[i] = pubkey$add(x[i], y[i])
# print(privkey$decrypt(c(x[i], y[i], res[i])))
res[i+1] = pubkey$add(x[i+1], y[i+1])
}
res
}
############################################################################
# We next define subtraction of two encrypted numbers: (enc(a1), enc(b1), c1) and (enc(a2), enc(b2), c2).
# We simply multiply the second encrypted number by -1 and then add it to the first encrypted number.
############
subenc = function(x, y) {
for (i in seq(1, length(x), 3)) {
# y[i:i+1] = pubkey$mult(y[i:i+1], -1) # this doesn't work, why?
y[i] = pubkey$mult(y[i], -1)
y[i+1] = pubkey$mult(y[i+1], -1)
}
addenc(x, y)
}
#############################################################################################
# We next define the multiplication of an encrypted number en = (enc(a),enc(b),c) by a scalar y
# Again, Paillier only allows y to be an integer. To deal with floating point numbers, we use
# en * y = en * (y_int + y_frac)
# = en * y_int + (enc(a) * floor(y_frac * 10^n), enc(b) * floor(y_frac * 10^n), n)
# where n is a parameter with larger n leading to better approximation. We use n = 3 below.
#######
smultenc = function(x,y) {
yint = floor(y * 1000)
res = x
for (i in seq(1, length(x), 3)) {
res[i] = pubkey$mult(x[i], yint)
res[i+1] = pubkey$mult(x[i+1], yint)
res[i+2] = x[i+2] + 3
}
res
}
################
# Here are some test cases
####
# five = encrypt(5.4)
# ten = encrypt(10.2)
# five
# ten
# decrypt(addenc(five,ten))
# decrypt(subenc(ten,five))
# decrypt(subenc(five,ten))
# decrypt(addenc(c(five,ten), c(ten,ten)))
#
# five = encrypt(5.4,1)
# ten = encrypt(10.2,2)
# decrypt(addenc(five,ten))
# decrypt(subenc(ten,five))
# decrypt(addenc(c(five,ten), c(ten,ten)))
#
# decrypt(smultenc(five, 4))
# decrypt(smultenc(ten, 6))
# decrypt(smultenc(c(five,ten), 4))
# decrypt(smultenc(c(five,ten), -4))Here are the functions to setup the parties in the computation.
x1 = c() # these 3 vectors belong to RE1
x2 = c()
x3 = c()
re1setup = function() {
x1 <<- encrypt(x[,1])
x2 <<- encrypt(x[,2])
x3 <<- encrypt(x[,3])
}
x4 = c() # this belong to RE2
re2setup = function() {
x4 <<- encrypt(x[,4])
}
mastersetup = function() {
re1setup()
re2setup()
}The function master_grad below computes the gradient: the difference between predicted values (encrypted) by REs and the true values
Note the gradient is unencrypted. We could encrypt it if necessary.
master_grad <- function(x) {
xd = decrypt(x)
xd - y
# subenc(xd, y)
}Function re1pred computes the partial predictions based on RE1 data only. The values are encrypted.
re1pred <- function() {
# return (theta1 * x[,1] + theta2 * x[,2] + theta3 * x[,3])
return (addenc(smultenc(x1, theta1),
addenc(smultenc(x2, theta2),
smultenc(x3, theta3)))
)
}RE2 takes the partial prediction of RE1 and then add its contribution to the prediction. The final prediction values are encrypted.
re2pred <- function(z) {
# return (z + theta4 * x[,4])
return (addenc(z, smultenc(x4, theta4)))
}We need a variable to keep track of progress.
cadj = c(0,0,0,0) # a variable we want to print to keep track of progressRE1 model update function and implements the gradient descent update formula. It takes the gradient (unencrypted) from the master and then adjust its own theta. The whole computation is assumed to take place within RE1 and no encryption is required.
re1update <- function(z) {
theta1 <<- theta1 - (alpha / n) * sum(z * x[,1])
theta2 <<- theta2 - (alpha / n) * sum(z * x[,2])
theta3 <<- theta3 - (alpha / n) * sum(z * x[,3])
cadj[1] <<- sum(z * x[,1])
cadj[2] <<- sum(z * x[,2])
cadj[3] <<- sum(z * x[,3])
}RE2 model update function which implements the gradient descent update formula.
re2update <- function(z) {
theta4 <<- theta4 - (alpha / n) * sum(z * x[,4])
cadj[4] <<- sum(z * x[,4])
}After all that, we can now run the privacy-preserving linear regression algorithm. The algorithm can be sensitive to the learning rate.
pp_lm_sgd <- function(iter, rate) {
theta1 <<- 0
theta2 <<- 0
theta3 <<- 0
theta4 <<- 0
alpha <<- rate
mastersetup() # encrypt the data and set up communication
for (i in 1:iter) {
p1 = re1pred() # partial prediction
px = re2pred(p1) # full prediction
grad = master_grad(px) # compute gradient based on difference between true values and predicted values
re1update(grad) # update models independently
re2update(grad)
print(cadj)
}
print(c(theta1,theta2,theta3,theta4))
}Results
First the standard result, not private.
model0 = lm_sgd(100, 0.1)## [1] -850.00000 118.25073 -48.15952 255.53613
## [1] -765.00000 106.44608 -40.56685 230.69506
## [1] -688.50000 95.87435 -34.00341 208.27524
## [1] -619.65000 86.40182 -28.33898 188.04018
## [1] -557.68500 77.90971 -23.45922 169.77652
## [1] -501.91650 70.29246 -19.26379 153.29177
## [1] -451.72485 63.45623 -15.66464 138.41227
## [1] -406.55236 57.31756 -12.58456 124.98134
## [1] -365.897128 51.802222 -9.955889 112.857635
## [1] -329.307416 46.844136 -7.719341 101.913607
## [1] -296.376674 42.384486 -5.823009 92.034188
## [1] -266.739007 38.370887 -4.221473 83.115553
## [1] -240.065106 34.756662 -2.875014 75.064024
## [1] -216.058595 31.500205 -1.748921 67.795075
## [1] -194.4527359 28.5644164 -0.8128824 61.2324350
## [1] -175.00746228 25.91619578 -0.04044792 55.30728536
## [1] -157.5067161 23.5260007 0.5914422 49.9575257
## [1] -141.756044 21.367451 1.102875 45.127119
## [1] -127.580440 19.416979 1.511329 40.765499
## [1] -114.82240 17.65352 1.83200 36.82703
## [1] -103.340156 16.058238 2.078083 33.270540
## [1] -93.006141 14.614279 2.261018 30.058861
## [1] -83.705527 13.306556 2.390715 27.158457
## [1] -75.334974 12.121558 2.475744 24.539060
## [1] -67.801477 11.047176 2.523503 22.173354
## [1] -61.021329 10.072556 2.540369 20.036686
## [1] -54.919196 9.187959 2.531827 18.106804
## [1] -49.427276 8.384644 2.502583 16.363624
## [1] -44.484549 7.654762 2.456670 14.789019
## [1] -40.03609 6.99126 2.39753 13.36663
## [1] -36.032485 6.387795 2.328094 12.081678
## [1] -32.429236 5.838664 2.250848 10.920840
## [1] -29.186312 5.338731 2.167895 9.872079
## [1] -26.267681 4.883371 2.081003 8.924531
## [1] -23.640913 4.468418 1.991654 8.068388
## [1] -21.276822 4.090112 1.901081 7.294797
## [1] -19.149140 3.745065 1.810302 6.595763
## [1] -17.234226 3.430215 1.720152 5.964070
## [1] -15.510803 3.142798 1.631309 5.393203
## [1] -13.959723 2.880316 1.544314 4.877278
## [1] -12.563751 2.640510 1.459595 4.410983
## [1] -11.307375 2.421334 1.377482 3.989522
## [1] -10.176638 2.220940 1.298221 3.608564
## [1] -9.158974 2.037649 1.221988 3.264197
## [1] -8.243077 1.869943 1.148902 2.952890
## [1] -7.418769 1.716444 1.079031 2.671454
## [1] -6.676892 1.575901 1.012403 2.417006
## [1] -6.0092029 1.4471792 0.9490125 2.1869464
## [1] -5.4082826 1.3292482 0.8888269 1.9789246
## [1] -4.8674544 1.2211708 0.8317909 1.7908184
## [1] -4.3807089 1.1220949 0.7778318 1.6207107
## [1] -3.9426380 1.0312455 0.7268632 1.4668697
## [1] -3.5483742 0.9479173 0.6787881 1.3277312
## [1] -3.1935368 0.8714678 0.6335017 1.2018819
## [1] -2.8741831 0.8013118 0.5908941 1.0880452
## [1] -2.5867648 0.7369158 0.5508516 0.9850675
## [1] -2.3280883 0.6777935 0.5132589 0.8919069
## [1] -2.0952795 0.6235009 0.4780000 0.8076216
## [1] -1.8857515 0.5736331 0.4449597 0.7313608
## [1] -1.6971764 0.5278201 0.4140243 0.6623556
## [1] -1.5274588 0.4857242 0.3850822 0.5999113
## [1] -1.3747129 0.4470366 0.3580247 0.5434000
## [1] -1.2372416 0.4114750 0.3327464 0.4922543
## [1] -1.1135174 0.3787814 0.3091455 0.4459611
## [1] -1.0021657 0.3487196 0.2871238 0.4040570
## [1] -0.9019491 0.3210735 0.2665872 0.3661229
## [1] -0.8117542 0.2956452 0.2474459 0.3317800
## [1] -0.7305788 0.2722535 0.2296140 0.3006858
## [1] -0.6575209 0.2507324 0.2130098 0.2725309
## [1] -0.5917688 0.2309298 0.1975557 0.2470353
## [1] -0.5325919 0.2127064 0.1831780 0.2239459
## [1] -0.4793327 0.1959341 0.1698071 0.2030339
## [1] -0.4313995 0.1804958 0.1573774 0.1840924
## [1] -0.3882595 0.1662839 0.1458265 0.1669341
## [1] -0.3494336 0.1531997 0.1350962 0.1513900
## [1] -0.3144902 0.1411526 0.1251312 0.1373069
## [1] -0.2830412 0.1300593 0.1158799 0.1245464
## [1] -0.2547371 0.1198436 0.1072936 0.1129831
## [1] -0.22926336 0.11043512 0.09932677 0.10250396
## [1] -0.20633703 0.10176956 0.09193665 0.09300629
## [1] -0.18570333 0.09378764 0.08508322 0.08439743
## [1] -0.16713299 0.08643494 0.07872903 0.07659346
## [1] -0.15041969 0.07966142 0.07283905 0.06951846
## [1] -0.13537772 0.07342108 0.06738056 0.06310371
## [1] -0.12183995 0.06767163 0.06232300 0.05728705
## [1] -0.10965596 0.06237417 0.05763786 0.05201218
## [1] -0.09869036 0.05749292 0.05329854 0.04722818
## [1] -0.08882132 0.05299496 0.04928024 0.04288892
## [1] -0.07993919 0.04885001 0.04555986 0.03895267
## [1] -0.07194527 0.04503021 0.04211588 0.03538164
## [1] -0.06475075 0.04150991 0.03892827 0.03214160
## [1] -0.05827567 0.03826550 0.03597840 0.02920157
## [1] -0.05244810 0.03527528 0.03324894 0.02653350
## [1] -0.04720329 0.03251922 0.03072377 0.02411196
## [1] -0.04248296 0.02997891 0.02838790 0.02191395
## [1] -0.03823467 0.02763741 0.02622744 0.01991862
## [1] -0.03441120 0.02547910 0.02422945 0.01810708
## [1] -0.03097008 0.02348958 0.02238193 0.01646223
## [1] -0.02787307 0.02165563 0.02067374 0.01496856
## [1] -0.02508577 0.01996504 0.01909455 0.01361203
## [1] 16.9995485 -2.4880167 0.3689761 -5.2744076Next the privacy preserving version.
model1 = pp_lm_sgd(100, 0.1)## [1] -850.00000 118.25073 -48.15952 255.53613
## [1] -765.00000 106.42593 -40.57662 230.65067
## [1] -688.50000 95.84965 -34.01919 208.25028
## [1] -619.65000 86.36412 -28.35340 188.04255
## [1] -557.70000 77.90824 -23.46557 169.78133
## [1] -501.95000 70.28198 -19.26387 153.27482
## [1] -451.75000 63.43096 -15.68613 138.37614
## [1] -406.55000 57.28054 -12.57338 124.98490
## [1] -365.90000 51.79797 -9.95896 112.85871
## [1] -329.300000 46.820055 -7.713993 101.901156
## [1] -296.40000 42.39851 -5.84545 92.01288
## [1] -266.75000 38.37015 -4.22447 83.09747
## [1] -240.100000 34.747150 -2.897562 75.058168
## [1] -216.050000 31.513309 -1.764233 67.794483
## [1] -194.450000 28.573349 -0.803082 61.211135
## [1] -175.00000000 25.88772662 -0.05365162 55.31072780
## [1] -157.5000000 23.5396075 0.5668839 49.9410515
## [1] -141.750000 21.381992 1.086892 45.106184
## [1] -127.600000 19.416239 1.507618 40.757125
## [1] -114.850000 17.643709 1.830308 36.844875
## [1] -103.350000 16.067119 2.057454 33.271433
## [1] -93.000000 14.628014 2.247512 30.036914
## [1] -83.700000 13.277394 2.409937 27.142676
## [1] -75.350000 12.124074 2.478064 24.538247
## [1] -67.80000 11.06196 2.52105 22.17610
## [1] -61.05000 10.04205 2.54835 20.05759
## [1] -54.900000 9.165055 2.543546 18.082008
## [1] -49.450000 8.379268 2.513601 16.348704
## [1] -44.500000 7.637043 2.469217 14.810042
## [1] -40.050000 6.990101 2.403430 13.366661
## [1] -36.050000 6.388080 2.324449 12.068920
## [1] -32.450000 5.830982 2.232274 10.916820
## [1] -29.200000 5.359709 2.167697 9.858756
## [1] -26.250000 4.893814 2.050382 8.948938
## [1] -23.650000 4.466104 1.971365 8.085514
## [1] -21.300000 4.086035 1.881647 7.269731
## [1] -19.150000 3.741432 1.827735 6.598343
## [1] -17.250000 3.404926 1.723577 5.977201
## [1] -15.500000 3.165060 1.599262 5.402340
## [1] -13.950000 2.865379 1.532157 4.876593
## [1] -12.550000 2.662339 1.444894 4.397127
## [1] -11.300000 2.408940 1.365842 3.968020
## [1] -10.200000 2.192367 1.323842 3.634309
## [1] -9.150000 2.035608 1.224632 3.251483
## [1] -8.250000 1.866675 1.171930 2.965412
## [1] -7.400000 1.707198 1.070228 2.680587
## [1] -6.700000 1.587265 1.008071 2.393157
## [1] -6.0000000 1.4646139 0.9434214 2.2037264
## [1] -5.400000 1.343322 0.880018 1.965296
## [1] -4.8500000 1.2206703 0.8153687 1.7758657
## [1] -4.4000000 1.1362039 0.7890178 1.6328301
## [1] -3.9500000 1.0121932 0.7231227 1.4923998
## [1] -3.5500000 0.9277268 0.6967718 1.3493642
## [1] -3.2000000 0.8527161 0.6214209 1.2075746
## [1] -2.8500000 0.8172497 0.5856143 1.0631798
## [1] -2.6000000 0.7314240 0.5580175 0.9691442
## [1] -2.3500000 0.6945983 0.5209650 0.8737494
## [1] -2.1000000 0.6074133 0.4921223 0.8287139
## [1] -1.9000000 0.5705876 0.4550698 0.7333191
## [1] -1.7000000 0.5337619 0.4180173 0.6379242
## [1] -1.5000000 0.4955769 0.3797189 0.5915294
## [1] -1.3500000 0.4573920 0.3414206 0.5451346
## [1] -1.2500000 0.4097512 0.3521222 0.4974939
## [1] -1.1000000 0.3715663 0.3138238 0.4510990
## [1] -1.0000000 0.3333813 0.2755254 0.4047042
## [1] -0.9000000 0.3347406 0.2767713 0.3557042
## [1] -0.8000000 0.2951963 0.2372270 0.3583094
## [1] -0.7500000 0.2475556 0.2479286 0.3106687
## [1] -0.6500000 0.2583706 0.2001745 0.2629146
## [1] -0.6000000 0.2093706 0.2096302 0.2642739
## [1] -0.5500000 0.2201856 0.1618761 0.2165197
## [1] -0.5000000 0.2201856 0.1618761 0.2165197
## [1] -0.4500000 0.1725449 0.1725777 0.1688790
## [1] -0.4000000 0.1820006 0.1235777 0.1701249
## [1] -0.3500000 0.1330006 0.1330335 0.1714842
## [1] -0.3000000 0.1343599 0.1342793 0.1224842
## [1] -0.3000000 0.1343599 0.1342793 0.1224842
## [1] -0.25000000 0.14381561 0.08527935 0.12373010
## [1] -0.20000000 0.09617490 0.09598096 0.07608939
## [1] -0.20000000 0.09617490 0.09598096 0.07608939
## [1] -0.20000000 0.09617490 0.09598096 0.07608939
## [1] -0.15000000 0.09617490 0.09598096 0.07608939
## [1] -0.15000000 0.10563062 0.04698096 0.07733527
## [1] -0.15000000 0.05663062 0.05643668 0.07869456
## [1] -0.10000000 0.05663062 0.05643668 0.07869456
## [1] -0.10000000 0.05798991 0.05768257 0.02969456
## [1] -0.10000000 0.05798991 0.05768257 0.02969456
## [1] -0.10000000 0.05798991 0.05768257 0.02969456
## [1] -0.10000000 0.05798991 0.05768257 0.02969456
## [1] -0.05000000 0.05798991 0.05768257 0.02969456
## [1] -0.050000000 0.067445632 0.008682572 0.030940450
## [1] -0.05000000 0.01844563 0.01813830 0.03229974
## [1] -0.05000000 0.01844563 0.01813830 0.03229974
## [1] -0.05000000 0.01844563 0.01813830 0.03229974
## [1] -0.05000000 0.01844563 0.01813830 0.03229974
## [1] -0.05000000 0.01844563 0.01813830 0.03229974
## [1] -0.05000000 0.01844563 0.01813830 0.03229974
## [1] -0.05000000 0.01844563 0.01813830 0.03229974
## [1] 1.101341e-13 1.844563e-02 1.813830e-02 3.229974e-02
## [1] 1.101341e-13 1.844563e-02 1.813830e-02 3.229974e-02
## [1] 17.0001000 -2.4874175 0.3696221 -5.2740215