lab:projects:01bayesian_differential_equations:using_winbugs_for_data_analysis
Table of Contents
WinBUGS
Generate data
## Generate data from a differential equation ## y' = .1y t<-1:100 y<-exp(.05*t)+rnorm(100) plot(t, y, type='l') dput(y) dput(seq(1,100,1))
library(deSolve)
N=100
##########################################################
## Simulate data from predator-prey model
## x'=4x-0.8xy
## y'=-3y+0.5xy
## sample size=N
##########################################################
##----------------------------------##
## differential function
##----------------------------------##
prey <- function(t, state, parameters){
with(as.list(c(state, parameters)),{
dx = r1*x + a1*x*y
dy = r2*y + a2*x*y
list(c(dx,dy))
}
)
}
##----------------------------------##
## Simulate data
##----------------------------------##
parameters <- c(r1 = 4, a1 = -0.8, r2 = -3, a2 = 0.5)
state <- c(y = 1, x = 1)
times <- seq(0,by=1,length=N)
out <- as.data.frame(ode(y=state,times=times,func=prey,parms=parameters))
t <- out$times
x <- out$x
y <- out$y
dput(x)
dput(y)
dput(seq(1,100,1))
##########################################################
## Simulate data from Lorenz System model
## x'=10(y-x)
## y'=28x-y-xz
## z'=xy-2.667z
##########################################################
##----------------------------------##
## differential function
##----------------------------------##
lorenz <- function(t, state, parameters){
with(as.list(c(state, parameters)),{
dx = sigma*(y-x)
dy = rho*x-y-x*z
dz = x*y-beta*z
list(c(dx,dy,dz))
}
)
}
##----------------------------------##
## Simulate data
##----------------------------------##
parameters <- c(sigma = 10, rho = 28, beta = 2.667)
state <- c(x = 1, y = 1, z = 1)
times <- seq(0,by=1,length=N)
out <- as.data.frame(ode(y=state,times=times,func=lorenz,parms=parameters))
t <- out$times
x <- out$x
y <- out$y
z <- out$z
dput(x)
dput(y)
dput(z)
dput(seq(1,100,1))
##########################################################
## Simulate data from second-order DE y''=ay'+ by
## x'=-0.05x-0.3y
## y'=x
##########################################################
##----------------------------------##
## differential function
##----------------------------------##
sde <- function(t, state, parameters){
with(as.list(c(state, parameters)),{
dy = x
dx = a*x + b*y
list(c(dy,dx))
}
)
}
##----------------------------------##
## Simulate data
##----------------------------------##
parameters <- c(a = -.05, b = -.3)
state <- c(y = 1, x = 1)
times <- seq(0,by=1,length=N)
out <- as.data.frame(ode(y=state,times=times,func=sde,parms=parameters))
t <- out$times
x <- out$x
y <- out$y
dput(x)
dput(y)
dput(seq(1,100,1))
A working example for using WinBUGS
## model for estimation
model{
x[1:n]<-ode(init, grid[1:n], D(C, t), origin, tol)
D(C, t)<- beta *C
init~dnorm(0, 1.0E-6)
beta~dnorm(0, .001)
for (i in 1:n){
y[i]~dnorm(x[i], pre_sig)
}
pre_sig~dgamma(.001,.001)
sig<-1/pre_sig
}
## starting values
list(init=1, beta=.01, pre_sig=1)
list(origin = 0,tol = 1.0E-5, n=100,
y=c(1.91628041660203, -0.524273077742236, -0.215600435761198, -0.499805150233081,
1.49453111255607, 0.583655583988035, 2.00575990676020, 0.18875828177211,
1.74052521837843, 1.68393418973157, 1.07077486358785, 2.58972948516908,
3.16935915557293, 2.97103951269921, 2.95602403187493, 3.63432922028808,
2.91978767331905, 3.11569655315135, 3.39107447913997, 2.08055461616841,
2.34575694771166, 3.47565554333372, 3.28952279896791, 4.67583057471815,
1.94275655180682, 4.06835083026542, 4.03351182488062, 2.78200223879326,
2.83930431068625, 5.88712054894954, 4.36022502134429, 5.40517202664213,
4.74178342482887, 5.96421070340728, 5.52087543323209, 5.11422707527591,
4.74625691650396, 6.33031808471635, 6.85447791362293, 6.78149810189367,
6.65735568069162, 7.39308294675777, 7.59117819110668, 8.4203629334473,
7.46277570811714, 9.66033841406203, 9.38958283255697, 11.8116039024227,
9.31989604165274, 13.9300422110149, 13.9180644766879, 14.2070208398851,
14.1666532236633, 14.0069727194287, 17.1637319051880, 15.1208449878726,
17.8597615695121, 19.0465188775196, 18.6270565772774, 19.9377454780847,
21.3873978722973, 22.0319381266234, 22.9802767837657, 23.0350057651524,
25.3284634957884, 26.9008196424777, 28.0335358743345, 30.3354821267959,
32.2451433186996, 33.3709985948714, 34.8909412106869, 36.5532679256467,
38.2340611159950, 39.7393723977161, 42.7416948441582, 46.4753853607269,
47.3762153340991, 48.8565016198419, 50.4596263427143, 56.6833600951975,
57.3587023289551, 60.9654050208031, 62.0093466124524, 66.8845540221531,
70.682644852117, 75.1740929847133, 76.0642985964311, 81.0932668873101,
85.3454738213016, 89.2134032374966, 95.1338072454866, 100.361827110599,
104.301504458540, 109.320129944412, 115.228561145994, 122.031062036943,
128.971236291414, 134.345773559877, 140.632205954776, 148.074984897580
),
grid=c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49,
50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65,
66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81,
82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97,
98, 99, 100)
)
## Predator-Prey model
model{
solution[1:n, 1:dim]<-ode(init[1:dim], grid[1:n], D(C[1:dim],t), origin, tol)
D(C[1], t)<-r1*C[1]-a1*C[1]*C[2]
D(C[2], t)<-(-r2*C[2])+a2*C[1]*C[2]
init[1]~dnorm(0, 1.0E-6)
init[2]~dnorm(0, 1.0E-6)
r1~dnorm(0, .001)
a1~dnorm(0, .001)
r2~dnorm(0, .001)
a2~dnorm(0, .001)
for(i in 1:n){
x[i]~dnorm(solution[i,1], pre_sig1)
y[i]~dnorm(solution[i,2], pre_sig2)
}
pre_sig1~dgamma(.001,.001)
pre_sig2~dgamma(.001,.001)
sig1<-1/pre_sig1
sig2<-1/pre_sig2
}
## starting values
list(init=c(1,1), r1=4.5, r2=3, a1=0.8, a2=0.1 pre_sig1=1, pre_sig2=1)
list(origin = 0, tol = 1.0E-5, n = 100, dim = 2,
x=c(1, 27.2493730144343, 0.367279274917532, 8.80300260698433, 0.378628361507895,
2.19531314365834, 6.28669346240878, 0.632871170911329, 19.6635431447406,
0.306877446623149, 5.06923837762363, 0.727917440566333, 1.30392265887201,
23.9865822593632, 0.43143441015614, 11.7836534609542, 0.322409034684567,
2.92979119689851, 2.54620350065872, 0.802970470981476, 24.5584831276482,
0.330441314319186, 6.81813568236758, 0.480674383174515, 1.71876757814502,
13.5212206527526, 0.524074521019393, 15.6450686251951, 0.303319695294992,
3.92951373159578, 1.19753929994834, 1.03647861238963, 27.318159227581,
0.374808756576609, 9.16920461239181, 0.36790122252684, 2.28532990298458,
5.51376445504854, 0.65389605721012, 20.3653277880607, 0.309128881483225,
5.28880962579076, 0.6792147613296, 1.35638065361539, 22.7399215187337,
0.443132183898499, 12.2907686019695, 0.317864342551642, 3.05906116100026,
2.25194054128866, 0.833466871187378, 25.2061603465802, 0.335873947998433,
7.13077207234425, 0.457774413009327, 1.79437532393675, 11.9369867155523,
0.541504194366999, 16.3223926289681, 0.302899442193289, 4.1150417494468,
1.08317476406158, 1.08051484149931, 27.2172145579245, 0.384047284572539,
9.60955716256897, 0.356925283431677, 2.39402595554884, 4.74303440037069,
0.679332127000742, 21.1805098548039, 0.312224044526731, 5.55368507718278,
0.629875072053153, 1.41973185570748, 21.1129811329163, 0.457342028303078,
12.8976755814484, 0.313495174230473, 3.21486076558199, 1.9650014453801,
0.870224058300198, 25.8803534972108, 0.34265414121663, 7.50608832632895,
0.434620578330906, 1.88527090292459, 10.2705870223502, 0.562496995702743,
17.1204409995128, 0.303075458774406, 4.33748339404646, 0.971542520551176,
1.13329046132159, 26.8390316883056, 0.395256872767998, 10.1337684590833,
0.346065461517574, 2.52380381741617, 4.00760093678675),
y=c(1, 4.43127474862561, 2.81627527837235, 0.466444972483609, 8.36851199691762,
0.593588570898793, 19.4215330063024, 1.48264431102916, 1.13402246843081,
4.32559649553310, 0.440163416520315, 12.6778987442410, 0.821490845479253,
10.1217663865694, 2.24430352588426, 0.55366147608418, 6.65444109465153,
0.516010241105529, 17.9853303527598, 1.19652676067639, 2.18343021051778,
3.43477030883979, 0.438442305857428, 10.1769541305800, 0.684365560613417,
17.1485251754162, 1.79316819626249, 0.755291257712084, 5.2805751585465,
0.464386869361937, 15.1365914777555, 0.972956331965342, 4.98956224776079,
2.72931569113157, 0.474501576618742, 8.10729118699338, 0.581228878452443,
19.4074375455635, 1.43815140383442, 1.22949433358674, 4.18619296605838,
0.438222330785728, 12.2927431872431, 0.799540391320019, 11.2094629911017,
2.17210121584380, 0.574080205329927, 6.43328927941991, 0.506968069962872,
17.5981068812464, 1.15985882406194, 2.44958203203251, 3.318360020774,
0.441243341547497, 9.8361302405205, 0.666805159303815, 17.8778292808577,
1.73324169683032, 0.80445994027665, 5.09464231730815, 0.458765573769936,
14.6750773340703, 0.942922720390777, 5.72083687580081, 2.63249982352102,
0.485232456570942, 7.81558443762768, 0.567659664823688, 19.2932244320959,
1.38872314843596, 1.35633666334647, 4.03114518096661, 0.436694275020096,
11.8596769999823, 0.775277314344897, 12.4850043800748, 2.09204084649579,
0.600744648923205, 6.187913737843, 0.497240463608204, 17.1296528841351,
1.11940181484789, 2.80778263892303, 3.18986812244019, 0.445521552326331,
9.45845878784726, 0.647624632800273, 18.5339665889217, 1.66741472464012,
0.869041843108993, 4.8904244595474, 0.453080028455255, 14.1548235571282,
0.910133018610582, 6.67239868568556, 2.52672006450802, 0.499436365165659,
7.4965948576234, 0.553106710980875, 19.0582491293477),
grid=c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49,
50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65,
66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81,
82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97,
98, 99, 100)
)
## Lorenz model
model{
solution[1:n, 1:dim]<-ode(init[1:dim], grid[1:n], D(C[1:dim],t), origin, tol)
D(C[1], t)<-sigma*(C[2]-C[1])
D(C[2], t)<-rho*C[1]-C[2]-C[1]*C[3]
D(C[3], t)<-C[1]*C[2]-beta*C[3]
init[1]~dnorm(0, 1.0E-6)
init[2]~dnorm(0, 1.0E-6)
init[3]~dnorm(0, 1.0E-6)
sigma~dnorm(0, .001)
rho~dnorm(0, .001)
beta~dnorm(0, .001)
for(i in 1:n){
x[i]~dnorm(solution[i,1], pre_sig1)
y[i]~dnorm(solution[i,2], pre_sig2)
z[i]~dnorm(solution[i,3], pre_sig3)
}
pre_sig1~dgamma(.001,.001)
pre_sig2~dgamma(.001,.001)
pre_sig3~dgamma(.001,.001)
sig1<-1/pre_sig1
sig2<-1/pre_sig2
sig3<-1/pre_sig3
}
## starting values
list(init=c(1,1,1), sigma=10, rho=25, beta=2, pre_sig1=1, pre_sig2=1, pre_sig3=1)
list(origin = 0, tol = 1.0E-5, n = 100, dim = 3,
x=c(1, -9.37913726491717, -8.1751938776253, -7.45437142188846,
-10.0939326783036, -6.51241922380599, -9.7377171857395, -7.67099854845567,
-7.0040180958596, -9.99321310551298, -4.8971679335613, -11.6762686830654,
-3.33655761920973, -5.73645909803875, 8.96227136468952, -4.31154401902471,
-6.4746616626519, -15.8117079676247, 1.14327999967735, -3.28243279946662,
-12.9795464243791, -2.36890023621321, -2.02912017570026, 13.3156473701275,
1.13311528541584, -4.08334051282310, -13.9587145999381, -1.24678348477577,
-1.54559968618750, 11.7378128801797, 16.0032208996031, -3.66986695102125,
-12.5643543219622, -3.28096457617624, -1.25221731447857, -11.0781552426542,
-3.52734813076968, 0.105123408903435, -6.64653177835597, -7.5643020165231,
-6.78678756642795, -8.1046176464895, -0.744743314821774, 7.60594822077389,
8.9946120628748, 2.10708586368598, -7.39598743384607, -8.36348611035768,
0.934184381877345, 1.02318472809510, -12.7682338219043, -2.69089512597560,
-4.73738680695567, 14.3718080196719, -1.81651688352748, 2.50849842230557,
-12.8877354505077, 6.45705184837316, 0.00915367982581274, -5.48362926477406,
-9.57973428746041, -3.79257025637243, -14.3403430424708, 14.6212215569764,
-3.60032146534769, -13.9224246905381, -15.4572437647567, 3.58755298215106,
13.8411849343665, 10.0980320271030, -2.66842499743106, -3.4817101572167,
9.01085980241026, 2.13993764211872, -3.02246620559846, 10.1646198056224,
1.70136151546515, -4.65563026631379, 11.6644655639484, 1.94761054497078,
-5.51637087066344, -15.6952757755013, 1.38054806925278, 3.11323332728716,
13.5078178203216, 1.88966531155075, -2.99688267003906, -13.3703153168628,
-2.01740214506712, 3.11151570344386, 13.7220002829346, 4.37634991952408,
0.0470067110156934, -5.76576057023238, -8.9757695087526, -4.94006986326112,
-12.9840659872793, -0.458352293683489, -1.4923698709404, 13.1226327224450
),
y=c(1, -8.35581545175765, -9.56532694829006, -6.18945406932969,
-10.9306722466254, -6.97814803699841, -7.69681385808116, -9.8522205879396,
-4.5819752753716, -12.9399795764831, -3.74688778964786, -14.8169043154570,
-1.98353436033592, -9.83444679617344, 16.5716401836511, -2.03571426117949,
-10.6739605609039, -14.4184575385487, -0.391585985775191, -5.40893978113351,
-12.9712962139661, -1.84268876206979, -3.66652433578305, 10.9189148516828,
-0.0726996578805162, -6.32923254286013, -11.9117629149880, 2.10222084323678,
-2.66186081028832, 17.7134906913287, 22.4133447462612, -3.77222876103467,
-16.7590498433815, 1.23895996451924, -4.80996876736793, -15.3589362796038,
-0.595220878117489, 0.261744229165418, -9.4359288432197, -3.62123134525295,
-10.2293837217872, -2.06152344629817, -0.802959072271164, 12.0736910549176,
1.27614157901210, -4.726286768814, -11.5394726329885, -1.27141063543053,
2.56008288515753, 1.84700112719347, -13.5861818040817, -2.15798655622168,
-8.41980464649944, 9.75428512134016, 3.27118912866921, 4.18325617251020,
-19.2201614037564, -1.93518896519080, 0.352945007334987, -8.01342133568397,
-4.69276526724961, -5.54136833931971, -18.0488650446110, 5.62881782652179,
-4.90538519999496, -17.6779655494433, -7.81630640873018, 4.77601705949645,
17.4624159590797, -0.643682663517818, -3.83111507032489, -6.59722893981189,
3.3726783703015, 3.00945926804811, -5.68762427642052, 4.52587287196414,
2.11568894212966, -8.44713140978458, 4.64911791142862, -4.08582182838139,
-9.01157048921513, -16.5778339310918, -1.04645921071250, 5.0422345898722,
15.1702342359289, -0.537622525744513, -4.94858858427873, -15.4300146003076,
0.434742726541585, 4.95243727901168, 16.8788840898063, -2.52329900875594,
0.287370857325400, -8.31324418154128, -4.40503860029847, -7.37223959660941,
-8.15819643010487, 2.07692933878203, -2.60892582954787, 17.2159906286838
),
z=c(1, 29.3643927525484, 24.6199493006891, 27.4376447301346, 27.9797573020866,
23.9176878567253, 30.8861597875724, 22.5915666699715, 28.3710135463757,
24.869966142624, 24.6717685000157, 27.2114685560438, 23.5755751132713,
13.8860505253708, 12.7492503398810, 25.7820441466289, 15.669099039319,
38.3541303837000, 22.1849935786710, 21.7105913715242, 32.9096993322520,
21.0400577731127, 10.4464298661644, 35.7220097965273, 21.4289625502751,
17.9062457311836, 36.2726215636924, 25.4408114594833, 13.8476372946582,
22.9717951548596, 30.3533552678601, 21.1446960991787, 27.185178815595,
27.9107444463712, 25.2940515270683, 24.550955048971, 26.0924084553989,
13.1485100923288, 20.7578938110669, 30.5579919172470, 18.4557541692614,
33.0323844336029, 16.3247751501462, 17.8343070997790, 35.250292460831,
29.8477731525447, 17.9994638032206, 34.1558853307776, 20.9319534784527,
10.1967301948454, 31.7329941719784, 21.4006189663649, 12.0275247818403,
39.0426649527917, 27.9044537822999, 13.2221047833747, 24.5889024895005,
33.6020423850752, 16.3251300283512, 19.8620833176024, 33.4789235191451,
16.9498661365248, 30.7679062672392, 42.5841023652175, 17.9763455737370,
30.0185664056687, 42.7800769119416, 18.3031165378794, 30.0522284896145,
38.5011557385177, 16.4199561187733, 8.76995480344762, 33.5577623021345,
16.2780368801897, 9.01461517783606, 34.7607323460141, 17.1657717802896,
11.2108438113767, 37.5364650763762, 28.9878745615758, 15.3603654829082,
36.0346261230805, 24.038040823855, 20.4206811090178, 31.8542150746185,
24.3211702467393, 20.7600859054553, 31.1894004944121, 24.4360617058219,
19.1075489911540, 30.4310930910248, 30.949480888224, 15.1454901489185,
20.2915494093508, 32.5608943463519, 17.3280111328285, 37.3948019228782,
23.7953998187151, 12.7113355331188, 28.2677158457294),
grid=c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49,
50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65,
66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81,
82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97,
98, 99, 100)
)
## second-order DE model
model{
solution[1:n, 1:dim]<-ode(init[1:dim], grid[1:n], D(C[1:dim],t), origin, tol)
D(C[1], t)<-a*C[1]+b*C[2]
D(C[2], t)<-C[1]
init[1]~dnorm(0, 1.0E-6)
init[2]~dnorm(0, 1.0E-6)
a~dnorm(0, .001)
b~dnorm(0, .001)
for(i in 1:n){
x[i]~dnorm(solution[i,1], pre_sig1)
y[i]~dnorm(solution[i,2], pre_sig2)
}
pre_sig1~dgamma(.001,.001)
pre_sig2~dgamma(.001,.001)
sig1<-1/pre_sig1
sig2<-1/pre_sig2
}
## starting values
list(init=c(1,1), a=-0.1, b=-0.1, pre_sig1=1, pre_sig2=1)
list(origin = 0, tol = 1.0E-5, n = 100, dim = 2,
x=c(1, 0.531527947130369, -0.065784099138918, -0.615189491620755,
-0.962236767878465, -1.01775506188670, -0.780118159550971, -0.331440353235473,
0.189941730456045, 0.631689509886521, 0.87162147273815, 0.851108542055394,
0.588706488549915, 0.17109670089686, -0.274973363252488, -0.620816275601627,
-0.772623730256733, -0.696535968054192, -0.425381355758178, -0.0460562785563337,
0.327911983110668, 0.59006219291898, 0.671035780483072, 0.556559965705839,
0.28883596042776, -0.048258801831552, -0.355141006782378, -0.545706426567409,
-0.571244041341174, -0.43251429277073, -0.177119390982121, 0.116366129585949,
0.362329502928236, 0.492895378933266, 0.476430945783167, 0.324805492480869,
0.0878817579314938, -0.162566714716592, -0.354407234997003, -0.435750801754642,
-0.388772598045408, -0.233137732780422, -0.0185599421672952,
0.190849780995452, 0.335581884211518, 0.377486861061577, 0.309620523814412,
0.156704132204032, -0.0334747351622453, -0.20482497178673, -0.309365267882392,
-0.320520510938333, -0.239661581734035, -0.0943515899104367,
0.070797613647228, 0.207688194311922, 0.278632414977405, 0.266600413067330,
0.179072759554841, 0.0447100665025318, -0.0958595284540634, -0.202217195112752,
-0.245678950013537, -0.216908915930778, -0.127640408935192, -0.00629996164859407,
0.110920301794381, 0.190769331970693, 0.212282657171761, 0.172165957335011,
0.08487345341411, -0.0223829615138093, -0.118021418540326, -0.175314627822382,
-0.179782351615696, -0.132725539427237, -0.0500863955047196,
0.0428159699492826, 0.118968079431046, 0.157454300783269, 0.149129040203315,
0.0986516654589856, 0.0224830012048783, -0.0563868269834599,
-0.115317941894303, -0.138465796604438, -0.120969525597589, -0.0698044810534205,
-0.00121414428723824, 0.064377437338166, 0.108398082542553, 0.119337302724661,
0.0956863143366764, 0.0458818951170348, -0.0145868127300932,
-0.0679433541269527, -0.0993076907470011, -0.100802341939692,
-0.0734572515948998, -0.0264826696460647),
y=c(1, 1.78347967759123, 2.01977289126153, 1.66813661363563, 0.857598812238843,
-0.158150142253254, -1.07922786098899, -1.64739350463864, -1.71771523216548,
-1.29440545819235, -0.522343646946024, 0.36113082077566, 1.09845836225460,
1.48634757235669, 1.43114484807504, 0.97021438498091, 0.254884017741135,
-0.498298439433085, -1.07254309460392, -1.31269842706697, -1.16652406203481,
-0.694577328272656, -0.0474290653801323, 0.581693053451526, 1.01412852032897,
1.13606151867514, 0.92783877228427, 0.464984452452852, -0.107993749439637,
-0.62220814203651, -0.933777901408383, -0.963670322476185, -0.717092366641987,
-0.277896308259809, 0.219185832140460, 0.629473984833002, 0.840112540503456,
0.80072643934597, 0.534749810709419, 0.129138145403612, -0.293545035318985,
-0.6118417102934, -0.740008756539104, -0.650740395161553, -0.380116575938532,
-0.0142128559450846, 0.337902080179523, 0.576413387207532, 0.638794500374811,
0.515833089826469, 0.251661660733976, -0.0714453628426446, -0.3584053876505,
-0.529088458504962, -0.540455922313865, -0.39703355721113, -0.147300404208550,
0.132289746880926, 0.360491072705445, 0.474686339138094, 0.447846560511128,
0.294508481420239, 0.0646013640154647, -0.172529220225885, -0.348858307544445,
-0.417030123991148, -0.362865736321981, -0.207789224543031, -0.00097761090008755,
0.196026465089611, 0.327480664713595, 0.359068315119108, 0.286645061142024,
0.135950237796817, -0.0461931150566349, -0.206270946034844, -0.299675817527363,
-0.303002940328447, -0.219696395850209, -0.0777561500495835,
0.0794502862617172, 0.206315410317506, 0.268114442436139, 0.250384617236215,
0.162064391145915, 0.0318034583083448, -0.101180992087561, -0.198804849422907,
-0.234932899438021, -0.202254174159701, -0.113450287183986, 0.00339874982145734,
0.113578902808821, 0.18597143911505, 0.201760695967728, 0.159201280732016,
0.07328525208429, -0.0293539414143821, -0.118610687708744, -0.169665048211669
),
grid=c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49,
50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65,
66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81,
82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97,
98, 99, 100)
)