Sunday, October 19, 2014

Tuning Laplaces Demon II

I am continuing with my trying all algorithms of Laplaces Demon. It is actually quite a bit more work than I expected but I do find that some of the things get clearer. Now that I am close to the end of calculating this second batch I learned that there is loads of adaptive algorithms. The point of those adaptations is not so much getting the correct posterior distribution, but rather getting enough information so one can set up the other algorithms which can get the desired posterior. For example, in this post DRAM is the adaptive version of DRM which form such a pairing of algorithms.
Given all that it may be that I will redo this same exercise with a different estimation, but that is yet to be decided.

Adaptive-Mixture Metropolis

No specs




Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Algorithm = "AMM")

Acceptance Rate: 0.284
Algorithm: Adaptive-Mixture Metropolis
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
   beta[1]    beta[2]
2.73756468 0.00197592

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
         All Stationary
Dbar  45.095     44.425
pD   234.487      2.231
DIC  279.582     46.656
Initial Values:
[1] -10   0

Iterations: 10000
Log(Marginal Likelihood): NA
Minutes of run-time: 0.05
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 500
Recommended Burn-In of Un-thinned Samples: 5000
Recommended Thinning: 150
Specs: (NOT SHOWN HERE)
Status is displayed every 100 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 1000
Thinning: 10


Summary of All Samples
                Mean          SD        MCSE       ESS          LB      Median
beta[1]  -10.8694827  1.66603784 0.154095363  57.79995 -15.5489329 -10.2511972
beta[2]    0.2682103  0.04423248 0.004039057  58.80563   0.2003859   0.2543845
Deviance  45.0951406 21.65580992 1.178504393 534.75650  42.5189305  43.4676853
LP       -31.3536988 10.82813264 0.589266937 534.75770 -33.8231778 -30.5333506
                  UB
beta[1]   -8.3168119
beta[2]    0.3923192
Deviance  50.0480146
LP       -30.0738495


Summary of Stationary Samples
               Mean        SD        MCSE      ESS          LB      Median
beta[1]  -11.574823 2.1300652 0.205608788 285.1749 -16.3533092 -11.3357688
beta[2]    0.286758 0.0553065 0.005306451 287.9282   0.1924482   0.2804307
Deviance  44.425140 2.1124891 0.191862970 191.2321  42.4735763  43.8636196
LP       -31.027498 1.0677294 0.589266937 190.2838 -34.0072489 -30.7386292
                  UB
beta[1]   -7.9334688
beta[2]    0.4128574
Deviance  50.4694327
LP       -30.0463562

Affine-Invariant Ensemble Sampler

It seems to go somewhere, then gets stuck without an exit.


Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Iterations = 20000, Status = 2000, Thinning = 35, Algorithm = "AIES",
    Specs = list(Nc = 16, Z = NULL, beta = 1.1, CPUs = 1, Packages = NULL,
        Dyn.libs = NULL))

Acceptance Rate: 0.9773
Algorithm: Affine-Invariant Ensemble Sampler
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
     beta[1]      beta[2]
0.5252284175 0.0004811633

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
        All Stationary
Dbar 43.004     43.005
pD    0.053      0.000
DIC  43.057     43.005
Initial Values:
[1] -10   0

Iterations: 20000
Log(Marginal Likelihood): NA
Minutes of run-time: 0.8
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 399
Recommended Burn-In of Un-thinned Samples: 13965
Recommended Thinning: 27
Specs: (NOT SHOWN HERE)
Status is displayed every 2000 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 571
Thinning: 35


Summary of All Samples
                Mean         SD        MCSE       ESS          LB      Median
beta[1]  -10.2521485 0.72528513 0.153054828  9.623682 -12.7424753  -9.9662793
beta[2]    0.2513389 0.01927108 0.004080774 11.791793   0.2404582   0.2438793
Deviance  43.0041950 0.32410474 0.046924334 74.412690  42.5190044  43.0023753
LP       -30.3005775 0.16647750 0.024331033 74.198162 -30.8672783 -30.2965116
                  UB
beta[1]   -9.8180273
beta[2]    0.3153106
Deviance  44.0738314
LP       -30.0671736


Summary of Stationary Samples
                Mean           SD         MCSE      ESS          LB      Median
beta[1]   -9.9558233 0.0082421078 2.797169e-03 12.56157  -9.9743833  -9.9550952
beta[2]    0.2436365 0.0001725992 5.836733e-05 12.63574   0.2433173   0.2436223
Deviance  43.0047636 0.0021518709 7.092913e-04 12.84743  43.0002894  43.0048552
LP       -30.2976030 0.0009940593 2.433103e-02 12.86979 -30.2995034 -30.2976427
                  UB
beta[1]   -9.9405874
beta[2]    0.2440232
Deviance  43.0088727
LP       -30.2955405

Componentwise Hit-And-Run Metropolis

This never was able to get to the target.

Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Iterations = 40000, Status = 2000, Thinning = 30, Algorithm = "CHARM")

Acceptance Rate: 0.31229
Algorithm: Componentwise Hit-And-Run Metropolis
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
    beta[1]     beta[2]
3.580895236 0.002467357

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
        All Stationary
Dbar 44.445     45.021
pD    2.023      2.256
DIC  46.468     47.278
Initial Values:
[1] -10   0

Iterations: 40000
Log(Marginal Likelihood): NA
Minutes of run-time: 0.18
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 1064
Recommended Burn-In of Un-thinned Samples: 31920
Recommended Thinning: 31
Specs: (NOT SHOWN HERE)
Status is displayed every 2000 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 1333
Thinning: 30


Summary of All Samples
                Mean         SD       MCSE      ESS          LB      Median
beta[1]  -10.9964257 1.89283881 0.49785194 13.06079 -14.8229746 -10.9766992
beta[2]    0.2717979 0.04913034 0.01300998 11.03343   0.1856506   0.2705021
Deviance  44.4449406 2.01148697 0.18601589 82.06199  42.4984709  43.8291949
LP       -31.0303916 1.00254773 0.09222924 83.71010 -33.6460481 -30.7196890
                  UB
beta[1]   -7.6255135
beta[2]    0.3698866
Deviance  49.6364683
LP       -30.0586484


Summary of Stationary Samples
                Mean         SD       MCSE       ESS        LB     Median
beta[1]   -9.5579858 1.34107513 0.62509204  4.739982 -12.03957  -9.313436
beta[2]    0.2340237 0.03463118 0.01639968  4.804878   0.18134   0.227444
Deviance  45.0214688 2.12434347 0.28430825 16.656844  42.51149  44.655282
LP       -31.3029682 1.05482519 0.09222924 17.255007 -33.92636 -31.139132
                 UB
beta[1]   -7.398433
beta[2]    0.297938
Deviance  50.284536
LP       -30.061842

Delayed Rejection Adaptive Metropolis

This is an interesting algorithm. One can see during sampling the algorithm shifts from a faster to a slower sampling approach. The same shift in gears is seen in the plot. Notice that it recommends thinning 90. In fact I had  it to the point of proposing a thinning of 1000. Since the manual also states on using DRAM as final algorithm: 'DRAM may be used if diminishing adaptation occurs and adaptation ceases effectively'. Given these texts and effects, I tried a different problem, starting with wrong initial values. Indeed, it was able to get close to the true values in all such runs.



Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Thinning = 30, Algorithm = "DRAM")

Acceptance Rate: 0.5221
Algorithm: Delayed Rejection Adaptive Metropolis
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
     beta[1]      beta[2]
11.556472479  0.007722216

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
             All Stationary
Dbar     470.735     48.093
pD   1803475.978     35.962
DIC  1803946.712     84.055
Initial Values:
[1] -10   0

Iterations: 10000
Log(Marginal Likelihood): NA
Minutes of run-time: 0.2
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 165
Recommended Burn-In of Un-thinned Samples: 4950
Recommended Thinning: 270
Specs: (NOT SHOWN HERE)
Status is displayed every 100 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 333
Thinning: 30


Summary of All Samples
                 Mean           SD         MCSE       ESS           LB
beta[1]   -11.7526275    2.3119401   0.34580302  43.51566   -17.027894
beta[2]     0.2000943    0.4891327   0.02860336 130.61590    -1.487342
Deviance  470.7346820 1899.1977136 105.17220909 100.44119    42.511693
LP       -244.1848392  949.6008219  52.58640512 100.44148 -3481.777935
              Median           UB
beta[1]  -11.6929022   -7.8194808
beta[2]    0.2842366    0.4423707
Deviance  44.3755257 6945.9261923
LP       -31.0009124  -30.0634427


Summary of Stationary Samples
                Mean         SD       MCSE ESS         LB      Median
beta[1]  -11.7250338 2.48894626  0.2213773 168 -17.044268 -11.6851217
beta[2]    0.2921779 0.06547124  0.0053019 168   0.166793   0.2909701
Deviance  48.0932958 8.48081577  0.7474430 168  42.527075  45.2507580
LP       -32.8641423 4.24023747 52.5864051 168 -47.454476 -31.4673576
                  UB
beta[1]   -7.5769778
beta[2]    0.4250974
Deviance  77.3040995
LP       -30.0696373

Delayed Rejection Metropolis

This algorithm has the instruction to use the covariance matrix from for instance DRAM. So I pulled those and the summary of stationary samples as input.


Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = c(-11.72,
    0.29), Covar = covar, Algorithm = "DRM")

Acceptance Rate: 0.5659
Algorithm: Delayed Rejection Metropolis
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
     beta[1]      beta[2]
11.556472479  0.007722216

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
         All Stationary
Dbar  48.417     48.417
pD    59.001     59.001
DIC  107.419    107.419
Initial Values:
[1] -11.72   0.29

Iterations: 10000
Log(Marginal Likelihood): -38.65114
Minutes of run-time: 0.09
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 0
Recommended Burn-In of Un-thinned Samples: 0
Recommended Thinning: 10
Specs: (NOT SHOWN HERE)
Status is displayed every 100 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 1000
Thinning: 10


Summary of All Samples
                Mean          SD        MCSE      ESS        LB      Median
beta[1]  -11.6326715  2.89304045 0.111808577 891.8417 -18.12638 -11.5067562
beta[2]    0.2883743  0.07495814 0.002893592 894.3397   0.13874   0.2834856
Deviance  48.4174842 10.86289496 0.377770754 897.9490  42.52784  44.7049759
LP       -33.0262590  5.43029899 0.188915825 897.4927 -47.25058 -31.1763877
                 UB
beta[1]   -6.014343
beta[2]    0.452027
Deviance  76.915884
LP       -30.075104


Summary of Stationary Samples
                Mean          SD        MCSE      ESS        LB      Median
beta[1]  -11.6326715  2.89304045 0.111808577 891.8417 -18.12638 -11.5067562
beta[2]    0.2883743  0.07495814 0.002893592 894.3397   0.13874   0.2834856
Deviance  48.4174842 10.86289496 0.377770754 897.9490  42.52784  44.7049759
LP       -33.0262590  5.43029899 0.188915825 897.4927 -47.25058 -31.1763877
                 UB
beta[1]   -6.014343
beta[2]    0.452027
Deviance  76.915884
LP       -30.075104

Differential Evolution Markov Chain

Following LP, one can see this algorithm shift its step to step towards the target distribution. The same is visible in the samples.


Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Iterations = 70000, Status = 2000, Thinning = 36, Algorithm = "DEMC",
    Specs = list(Nc = 3, Z = NULL, gamma = 0, w = 0.1))

Acceptance Rate: 0.94571
Algorithm: Differential Evolution Markov Chain
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
    beta[1]     beta[2]
90.26633832  0.04206898

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
          All Stationary
Dbar   89.209     43.944
pD   5238.430      1.706
DIC  5327.639     45.650
Initial Values:
[1] -10   0

Iterations: 70000
Log(Marginal Likelihood): NA
Minutes of run-time: 0.73
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 1164
Recommended Burn-In of Un-thinned Samples: 41904
Recommended Thinning: 32
Specs: (NOT SHOWN HERE)
Status is displayed every 2000 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 1944
Thinning: 36


Summary of All Samples
               Mean          SD        MCSE      ESS           LB      Median
beta[1]  -17.482864   9.5017889  2.35451645 2.787539  -36.9570864 -13.2979369
beta[2]    0.410902   0.2049482  0.04994949 4.145223    0.1652045   0.3204944
Deviance  89.209031 102.3565307 23.04234800 7.261840   42.4986747  44.6939882
LP       -53.548197  51.3373427 11.56914508 7.198904 -167.9946428 -31.1456515
                  UB
beta[1]   -7.1204957
beta[2]    0.7804724
Deviance 317.1315784
LP       -30.0563707


Summary of Stationary Samples
                Mean        SD         MCSE       ESS          LB      Median
beta[1]  -11.9431454 1.7007792  0.215271093 125.13410 -15.7999987 -11.9807022
beta[2]    0.2969431 0.0441033  0.005730276 118.23767   0.2259635   0.2946702
Deviance  43.9443086 1.8471515  0.371329880  63.98536  42.4849394  43.3846382
LP       -30.7905955 0.9340381 11.569145079  63.51772 -33.2807813 -30.5163527
                  UB
beta[1]   -9.0646733
beta[2]    0.4059097
Deviance  48.8635918
LP       -30.0522792

Elliptical Slice Sampler

Manual states. 'This algorithm is applicable only to models in which the prior mean of all parameters is zero.' That is true for my prior, yet I am not impressed at all. Maybe I should be centering or such, but the current formulation was not a success
Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Iterations = 60000, Status = 2000, Thinning = 1000, Algorithm = "ESS")

Acceptance Rate: 1
Algorithm: Elliptical Slice Sampler
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
    beta[1]     beta[2]
1.514016386 0.001094917

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
        All Stationary
Dbar 53.903     53.806
pD   11.574     13.346
DIC  65.477     67.152
Initial Values:
[1] -10   0

Iterations: 60000
Log(Marginal Likelihood): NA
Minutes of run-time: 0.77
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 18
Recommended Burn-In of Un-thinned Samples: 18000
Recommended Thinning: 1000
Specs: (NOT SHOWN HERE)
Status is displayed every 2000 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 60
Thinning: 1000


Summary of All Samples
                Mean         SD        MCSE      ESS           LB     Median
beta[1]   -5.9519978 1.12538724 0.199063822 34.62487  -8.11739884  -5.904983
beta[2]    0.1419102 0.02788798 0.004825592 38.79318   0.09329411   0.141184
Deviance  53.9025233 4.81129661 0.854043909 46.96804  46.49740770  53.653605
LP       -35.7152403 2.39947768 0.425934607 46.93833 -41.22733361 -35.594423
                  UB
beta[1]   -3.9669525
beta[2]    0.1932619
Deviance  64.9487765
LP       -32.0237607


Summary of Stationary Samples
                Mean         SD        MCSE      ESS           LB      Median
beta[1]   -5.9962946 1.24391453 0.253903583 22.52123  -8.27636391  -5.9679392
beta[2]    0.1430514 0.03108438 0.006168722 27.21658   0.09456836   0.1467105
Deviance  53.8060933 5.16634523 1.088725764 34.31438  46.18528394  53.6113477
LP       -35.6674227 2.57618453 0.425934607 34.28162 -40.73404524 -35.5728340
                  UB
beta[1]   -4.0287456
beta[2]    0.1938871
Deviance  63.9614942
LP       -31.8687620

Gibbs Sampler

This needs derivatives, hence skipped.

Griddy Gibbs

This takes a grid from which a density is estimated and on which sampling is based. It may be a bit difficult for this grid, since the two parameters have different scales and the same grid is used. With only two parameters it was possible to take a rather high value for the number of grid points. Even so, I am not so happy with the final outcome.
Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Iterations = 30000, Status = 2000, Thinning = 100, Algorithm = "GG",
    Specs = list(Grid = seq(from = -0.25, to = 0.25, len = 13),
        dparm = NULL, CPUs = 1, Packages = NULL, Dyn.libs = NULL))

Acceptance Rate: 1
Algorithm: Griddy-Gibbs
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
     beta[1]      beta[2]
11.378198005  0.008486228

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
          All Stationary
Dbar   66.161     66.161
pD   1339.075   1339.075
DIC  1405.236   1405.236
Initial Values:
[1] -10   0

Iterations: 30000
Log(Marginal Likelihood): NA
Minutes of run-time: 2.09
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 0
Recommended Burn-In of Un-thinned Samples: 0
Recommended Thinning: 900
Specs: (NOT SHOWN HERE)
Status is displayed every 2000 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 300
Thinning: 100


Summary of All Samples
              Mean         SD       MCSE      ESS            LB      Median
beta[1]  -11.00845  3.3782928 0.77873105  23.0348  -18.26815566 -10.7755255
beta[2]    0.27170  0.0909315 0.01994284  30.0812    0.09175425   0.2612613
Deviance  66.16122 51.7508405 2.84979150 300.0000   42.82591844  50.8991200
LP       -41.89256 25.8754878 1.42498409 300.0000 -139.85980754 -34.2665415
                 UB
beta[1]   -4.870858
beta[2]    0.450951
Deviance 262.096671
LP       -30.229348


Summary of Stationary Samples
              Mean         SD       MCSE      ESS            LB      Median
beta[1]  -11.00845  3.3782928 0.77873105  23.0348  -18.26815566 -10.7755255
beta[2]    0.27170  0.0909315 0.01994284  30.0812    0.09175425   0.2612613
Deviance  66.16122 51.7508405 2.84979150 300.0000   42.82591844  50.8991200
LP       -41.89256 25.8754878 1.42498409 300.0000 -139.85980754 -34.2665415
                 UB
beta[1]   -4.870858
beta[2]    0.450951
Deviance 262.096671
LP       -30.229348


Hamiltonian Monte Carlo
A set was of specs was found. Acceptance rate is a bit high compared to the manual.

Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Thinning = 100, Algorithm = "HMC", Specs = list(epsilon = 0.9 *
        c(0.1, 0.01), L = 11))

Acceptance Rate: 0.8385
Algorithm: Hamiltonian Monte Carlo
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
    beta[1]     beta[2]
3.515108412 0.003083421

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
        All Stationary
Dbar 44.429     44.562
pD    1.941      1.869
DIC  46.369     46.431
Initial Values:
[1] -10   0

Iterations: 10000
Log(Marginal Likelihood): NA
Minutes of run-time: 0.59
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 80
Recommended Burn-In of Un-thinned Samples: 8000
Recommended Thinning: 100
Specs: (NOT SHOWN HERE)
Status is displayed every 100 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 100
Thinning: 100


Summary of All Samples
               Mean         SD        MCSE ESS          LB      Median
beta[1]  -11.421073 1.87894067 0.242559120 100 -15.4741232 -11.3104311
beta[2]    0.283175 0.04808956 0.006413119 100   0.1975276   0.2818055
Deviance  44.428764 1.97004886 0.159856183 100  42.5400966  43.7163906
LP       -31.027024 0.98807551 0.080511201 100 -33.5265769 -30.6632451
                  UB
beta[1]   -7.9608191
beta[2]    0.3807289
Deviance  49.3945952
LP       -30.0829425


Summary of Stationary Samples
                Mean         SD        MCSE ESS          LB      Median
beta[1]  -11.0974590 1.92886822 0.226325971  20 -15.4741232 -10.9792610
beta[2]    0.2750898 0.04775153 0.005911688  20   0.2034193   0.2713854
Deviance  44.5623988 1.93322645 0.408444645  20  42.5740058  44.0794655
LP       -31.0902147 0.97037005 0.080511201  20 -33.0194587 -30.8456818
                  UB
beta[1]   -8.2355095
beta[2]    0.3807289
Deviance  48.3972203
LP       -30.0962034

Another set of specs


Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Thinning = 100, Algorithm = "HMC", Specs = list(epsilon = 3 *
        c(0.1, 0.001), L = 18))

Acceptance Rate: 0.8855
Algorithm: Hamiltonian Monte Carlo
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
    beta[1]     beta[2]
3.640714435 0.003207219

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
        All Stationary
Dbar 44.404     44.404
pD    2.051      2.051
DIC  46.455     46.455
Initial Values:
[1] -10   0

Iterations: 10000
Log(Marginal Likelihood): NA
Minutes of run-time: 0.96
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 0
Recommended Burn-In of Un-thinned Samples: 0
Recommended Thinning: 100
Specs: (NOT SHOWN HERE)
Status is displayed every 100 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 100
Thinning: 100


Summary of All Samples
                Mean         SD        MCSE ESS          LB      Median
beta[1]  -11.5949171 1.91103354 0.200246790 100 -15.6570246 -11.5727273
beta[2]    0.2867121 0.04916803 0.005146306 100   0.2097083   0.2865395
Deviance  44.4043210 2.02528350 0.193624072 100  42.4813611  43.7159364
LP       -31.0168639 1.01912132 0.097084186 100 -33.7046786 -30.6665144
                  UB
beta[1]   -8.4758710
beta[2]    0.3936658
Deviance  49.8533556
LP       -30.0520014


Summary of Stationary Samples
                Mean         SD        MCSE ESS          LB      Median
beta[1]  -11.5949171 1.91103354 0.200246790 100 -15.6570246 -11.5727273
beta[2]    0.2867121 0.04916803 0.005146306 100   0.2097083   0.2865395
Deviance  44.4043210 2.02528350 0.193624072 100  42.4813611  43.7159364
LP       -31.0168639 1.01912132 0.097084186 100 -33.7046786 -30.6665144
                  UB
beta[1]   -8.4758710
beta[2]    0.3936658
Deviance  49.8533556
LP       -30.0520014

Sunday, October 12, 2014

Tuning LaplacesDemon

I was continuing with my Bayesian algorithms in R exercise. For these exercises I port SAS PROC MCMC examples to the various R solutions. However, the next example was logit model and that's just too simple, especially after last week's Jacobian for the Box-Cox transformation. Examples for logit and probit models abound on the web. Hence I took the opportunity to experiment a bit with the various algorithms which LaplacesDemon offers on the logit model. I did nine of them, so I might do some more next week.
One thing I noticed for all algorithms. LaplacesDemon has an in build calculation which chooses the length of the burn in and suggests a thinning. This does not work very well. Sometimes it seems for may samples on the target distribution, yet indicates all samples are burn in. Other times it thinks the burn in is done, yet does not find the expected results. Running several chains and human decisions to be more efficient than these long chains with unclear decisions.
It is a bit scary and disappointing that even though this is a very small and simple model and the number of samples is sometimes large, the posterior sometimes does not contain the frequentist solution. How would I know the answer is wrong for more complex models, where the answer is not so obvious?

Data

Data are from example 3 of SAS Proc MCMC. The PROC MCMC stimates are -11.77 and 0.292 respectively. It is a logistic model. The glm() estimates are thus:
Call:  glm(formula = cbind(y, n - y) ~ x, family = binomial, data = set2)

Coefficients:
(Intercept)            x 
   -11.2736       0.2793 

Degrees of Freedom: 19 Total (i.e. Null);  18 Residual
Null Deviance:        71.8
Residual Deviance: 15.25     AIC: 46.44

Setup

For each algorithm I started with the default number of samples and thinning. Based on the output the thinning and number of samples would be adapted until such point that the thinning was at least as large as the recommended thinning and there was a decent number of samples after burn in.

MWG


MWG needed 80000 samples before it gave answers. Notice that beta[2] was not mixing very well. The estimates seem off.


Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Iterations = 80000, Status = 2000, Thinning = 30)

Acceptance Rate: 0.24453
Algorithm: Metropolis-within-Gibbs
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
 beta[1]  beta[2]
2.835066 2.835066

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
        All Stationary
Dbar 45.242     45.360
pD    1.966      1.955
DIC  47.208     47.315
Initial Values:
[1] -10   0

Iterations: 80000
Log(Marginal Likelihood): -23.51154
Minutes of run-time: 0.35
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 1862
Recommended Burn-In of Un-thinned Samples: 55860
Recommended Thinning: 34
Specs: (NOT SHOWN HERE)
Status is displayed every 2000 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 2666
Thinning: 30


Summary of All Samples
                Mean         SD        MCSE       ESS          LB      Median
beta[1]   -9.0326306 1.04621897 0.220388172  7.322802 -11.2499597  -8.9042795
beta[2]    0.2209152 0.02658052 0.005719566  7.354199   0.1730201   0.2174406
Deviance  45.2421985 1.98273668 0.194161638 22.653492  42.5636767  44.7799650
LP       -31.4080976 0.98541260 0.095926332 23.196296 -33.8514305 -31.1737818
                  UB
beta[1]   -7.1874098
beta[2]    0.2772263
Deviance  50.1501559
LP       -30.0892713


Summary of Stationary Samples
                Mean         SD        MCSE      ESS          LB      Median
beta[1]   -8.9359418 0.99836162 0.349118267 3.269911 -10.6298830  -8.8780160
beta[2]    0.2188631 0.02541074 0.009077136 2.979785   0.1818618   0.2180453
Deviance  45.3601829 1.97741388 0.197527290 5.636751  42.7114169  45.2011226
LP       -31.4661714 0.98254393 0.095926332 5.760722 -33.7158111 -31.3894625
                  UB
beta[1]   -7.4135378
beta[2]    0.2609228
Deviance  49.8800864
LP       -30.1541170


HARM

No specs

I ended with 160000 samples. From the plot I'd say that is somewhat overkill.

Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Iterations = 160000, Status = 2000, Thinning = 36, Algorithm = "HARM")

Acceptance Rate: 0.06213
Algorithm: Hit-And-Run Metropolis
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
    beta[1]     beta[2]
2.863159336 0.001961247

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
        All Stationary
Dbar 44.449     44.449
pD    1.704      1.704
DIC  46.153     46.153
Initial Values:
[1] -10   0

Iterations: 160000
Log(Marginal Likelihood): -23.13795
Minutes of run-time: 0.38
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 0
Recommended Burn-In of Un-thinned Samples: 0
Recommended Thinning: 36
Specs: (NOT SHOWN HERE)
Status is displayed every 2000 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 4444
Thinning: 36


Summary of All Samples
                Mean         SD        MCSE       ESS          LB     Median
beta[1]  -12.5516619 1.69184474 0.251682993  43.36212 -15.7525684 -12.537176
beta[2]    0.3121025 0.04404281 0.006471107  46.72382   0.2306382   0.312301
Deviance  44.4489336 1.84619493 0.166569318 221.64136  42.5029914  43.927156
LP       -31.0503518 0.93299371 0.085590782 211.80050 -33.5115921 -30.783807
                  UB
beta[1]   -9.3734152
beta[2]    0.3968852
Deviance  49.3254872
LP       -30.0617562


Summary of Stationary Samples
                Mean         SD        MCSE       ESS          LB     Median
beta[1]  -12.5516619 1.69184474 0.251682993  43.36212 -15.7525684 -12.537176
beta[2]    0.3121025 0.04404281 0.006471107  46.72382   0.2306382   0.312301
Deviance  44.4489336 1.84619493 0.166569318 221.64136  42.5029914  43.927156
LP       -31.0503518 0.93299371 0.085590782 211.80050 -33.5115921 -30.783807
                  UB
beta[1]   -9.3734152
beta[2]    0.3968852
Deviance  49.3254872
LP       -30.0617562

Specs: list(alpha.star=0.234, B=NULL)

This sometimes got stuck, making all samples the same. Such faulty runs are repeated. In addition it seems the algorithm is not able to detect a good run. When I came to a whopping 1280000 samples I decided enough is enough. 



Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Iterations = 1280000, Status = 8000, Thinning = 42, Algorithm = "HARM",
    Specs = list(alpha.star = 0.234, B = NULL))

Acceptance Rate: 0.23416
Algorithm: Hit-And-Run Metropolis
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
    beta[1]     beta[2]
3.639417089 0.002441388

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
        All Stationary
Dbar 44.359     44.328
pD    1.597      1.532
DIC  45.956     45.860
Initial Values:
[1] -10   0

Iterations: 1280000
Log(Marginal Likelihood): -23.38909
Minutes of run-time: 2.91
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 6094
Recommended Burn-In of Un-thinned Samples: 255948
Recommended Thinning: 44
Specs: (NOT SHOWN HERE)
Status is displayed every 8000 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 30476
Thinning: 42


Summary of All Samples
                Mean         SD        MCSE       ESS         LB      Median
beta[1]  -11.8540232 1.90772737 0.141469541  28.83441 -15.618594 -11.8276283
beta[2]    0.2941959 0.04938246 0.003631577  29.75273   0.206555   0.2931946
Deviance  44.3593200 1.78715332 0.031531123 190.88005  42.500061  43.8415966
LP       -30.9974154 0.90059957 0.016062596 182.03381 -33.375122 -30.7354077
                 UB
beta[1]   -8.489146
beta[2]    0.392589
Deviance  49.076360
LP       -30.058998


Summary of Stationary Samples
                Mean         SD        MCSE       ESS          LB      Median
beta[1]  -11.7694571 1.87121305 0.153383946  24.00847 -15.4323928 -11.8230837
beta[2]    0.2920393 0.04844789 0.003937189  24.51596   0.2027361   0.2930413
Deviance  44.3276492 1.75044211 0.033218178 182.18810  42.5044563  43.8258628
LP       -30.9805115 0.87993039 0.016062596 176.04438 -33.2870801 -30.7284209
                  UB
beta[1]   -8.3411645
beta[2]    0.3869563
Deviance  48.9102683
LP       -30.0612554

Adaptive Directional Metropolis-within-Gibbs

no specs

This algorithm seemed to be able to get stuck outside the target region. As other algorithms it seems that increasing the thinning only leads to suggestions for further increase.
 

Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Iterations = 250000, Status = 2000, Thinning = 35, Algorithm = "ADMG")

Acceptance Rate: 0.0841
Algorithm: Adaptive Directional Metropolis-within-Gibbs
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
    beta[1]     beta[2]
3.694667553 0.002472681

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
        All Stationary
Dbar 44.362     44.362
pD    1.790      1.790
DIC  46.152     46.152
Initial Values:
[1] -10   0

Iterations: 250000
Log(Marginal Likelihood): NA
Minutes of run-time: 1.82
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 0
Recommended Burn-In of Un-thinned Samples: 0
Recommended Thinning: 38
Specs: (NOT SHOWN HERE)
Status is displayed every 2000 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 7142
Thinning: 35


Summary of All Samples
                Mean         SD        MCSE       ESS          LB      Median
beta[1]  -11.5514694 1.92227432 0.252091958  40.70347 -15.7178588 -11.4718511
beta[2]    0.2863682 0.04966254 0.006471524  40.98657   0.1965967   0.2829988
Deviance  44.3618406 1.89210234 0.068523057 227.09714  42.4858150  43.7976307
LP       -30.9951604 0.95027797 0.034622622 221.64249 -33.4547190 -30.7101949
                  UB
beta[1]   -8.0953413
beta[2]    0.3936748
Deviance  49.2894939
LP       -30.0514927


Summary of Stationary Samples
                Mean         SD        MCSE       ESS          LB      Median
beta[1]  -11.5514694 1.92227432 0.252091958  40.70347 -15.7178588 -11.4718511
beta[2]    0.2863682 0.04966254 0.006471524  40.98657   0.1965967   0.2829988
Deviance  44.3618406 1.89210234 0.068523057 227.09714  42.4858150  43.7976307
LP       -30.9951604 0.95027797 0.034622622 221.64249 -33.4547190 -30.7101949
                  UB
beta[1]   -8.0953413
beta[2]    0.3936748
Deviance  49.2894939
LP       -30.0514927

Specs: list(n = 0, Periodicity = 100)

This seems to be a run which does not mix very well in beta[1]. The posteriour is not on the target either. However, I do not see an obvious characteristic from which to notice the latter. 

Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Iterations = 64000, Status = 2000, Thinning = 33, Algorithm = "ADMG",
    Specs = list(n = 0, Periodicity = 100))

Acceptance Rate: 0.09229
Algorithm: Adaptive Directional Metropolis-within-Gibbs
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
    beta[1]     beta[2]
4.718765252 0.003150388

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
        All Stationary
Dbar 44.625     44.368
pD    1.848      0.914
DIC  46.473     45.282
Initial Values:
[1] -10   0

Iterations: 64000
Log(Marginal Likelihood): NA
Minutes of run-time: 0.36
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 1544
Recommended Burn-In of Un-thinned Samples: 50952
Recommended Thinning: 32
Specs: (NOT SHOWN HERE)
Status is displayed every 2000 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 1939
Thinning: 33


Summary of All Samples
               Mean         SD       MCSE      ESS          LB      Median
beta[1]  -11.314847 2.17303791 0.51002442 11.27064 -15.1041582 -11.4860852
beta[2]    0.279974 0.05601643 0.01310839 10.57305   0.1747337   0.2834938
Deviance  44.625141 1.92250413 0.19216964 59.12247  42.5136669  44.1268753
LP       -31.124617 0.95756642 0.09526603 60.52000 -33.5147916 -30.8800217
                  UB
beta[1]   -7.2420670
beta[2]    0.3795564
Deviance  49.4135552
LP       -30.0677263


Summary of Stationary Samples
                Mean         SD        MCSE       ESS          LB      Median
beta[1]  -13.3566117 0.72016481 0.299751317  5.473664 -15.0713103 -13.1769958
beta[2]    0.3325182 0.01971691 0.007968755  6.324196   0.3014903   0.3294776
Deviance  44.3682834 1.35175074 0.137999589 30.199161  42.7948298  44.0306538
LP       -31.0192877 0.68081646 0.095266030 28.411038 -32.6277457 -30.8439440
                  UB
beta[1]  -12.2282779
beta[2]    0.3766297
Deviance  47.5762853
LP       -30.2209601

Adaptive Griddy-Gibbs

A thinning of 1000 was proposed.

Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Iterations = 80000, Status = 2000, Thinning = 1000, Algorithm = "AGG",
    Specs = list(Grid = GaussHermiteQuadRule(3)$nodes, dparm = NULL,
        smax = Inf, CPUs = 1, Packages = NULL, Dyn.libs = NULL))

Acceptance Rate: 1
Algorithm: Adaptive Griddy-Gibbs
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
    beta[1]     beta[2]
19.42384441  0.01546567

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
          All Stationary
Dbar   53.739     53.739
pD   2320.547   2320.547
DIC  2374.286   2374.286
Initial Values:
[1] -10   0

Iterations: 80000
Log(Marginal Likelihood): NA
Minutes of run-time: 2.73
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 0
Recommended Burn-In of Un-thinned Samples: 0
Recommended Thinning: 1000
Specs: (NOT SHOWN HERE)
Status is displayed every 2000 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 80
Thinning: 1000


Summary of All Samples
               Mean         SD       MCSE ESS          LB      Median
beta[1]  -10.922145  4.4338567 0.50560810  80 -16.5329510 -11.2359881
beta[2]    0.269604  0.1214608 0.01367873  80   0.1295182   0.2783504
Deviance  53.738973 68.1255818 7.52694767  80  42.5397709  44.3774010
LP       -35.684516 34.0782560 3.76525972  80 -40.1716354 -31.0031726
                  UB
beta[1]   -5.6863543
beta[2]    0.4200815
Deviance  62.5858557
LP       -30.0729465


Summary of Stationary Samples
               Mean         SD       MCSE ESS          LB      Median
beta[1]  -10.922145  4.4338567 0.50560810  80 -16.5329510 -11.2359881
beta[2]    0.269604  0.1214608 0.01367873  80   0.1295182   0.2783504
Deviance  53.738973 68.1255818 7.52694767  80  42.5397709  44.3774010
LP       -35.684516 34.0782560 3.76525972  80 -40.1716354 -31.0031726
                  UB
beta[1]   -5.6863543
beta[2]    0.4200815
Deviance  62.5858557
LP       -30.0729465

Adaptive Hamiltonian Monte Carlo


Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Iterations = 150000, Status = 2000, Thinning = 36, Algorithm = "AHMC",
    Specs = list(epsilon = rep(0.02, length(Initial.Values)),
        L = 2, Periodicity = 10))

Acceptance Rate: 0.36918
Algorithm: Adaptive Hamiltonian Monte Carlo
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
     beta[1]      beta[2]
1.1442983280 0.0008202951

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
        All Stationary
Dbar 43.736     43.682
pD    1.214      1.272
DIC  44.950     44.954
Initial Values:
[1] -10   0

Iterations: 150000
Log(Marginal Likelihood): NA
Minutes of run-time: 1.74
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 2912
Recommended Burn-In of Un-thinned Samples: 104832
Recommended Thinning: 36
Specs: (NOT SHOWN HERE)
Status is displayed every 2000 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 4166
Thinning: 36


Summary of All Samples
                Mean         SD        MCSE        ESS          LB      Median
beta[1]  -11.4418192 1.06961398 0.191704524   7.357239 -14.1113189 -11.4442594
beta[2]    0.2837334 0.02830491 0.004949075   9.071151   0.2346106   0.2835194
Deviance  43.7356469 1.55816164 0.036118007 526.351697  42.4689017  43.1991684
LP       -30.6795260 0.78105728 0.018367981 492.409683 -32.7826226 -30.4108318
                  UB
beta[1]   -9.6350028
beta[2]    0.3515657
Deviance  47.9431739
LP       -30.0443227


Summary of Stationary Samples
                Mean         SD       MCSE        ESS          LB      Median
beta[1]  -11.9854873 0.56317415 0.15499286   8.169076 -12.9867511 -11.9369396
beta[2]    0.2982346 0.01576288 0.00396860  11.739663   0.2675003   0.2982497
Deviance  43.6824281 1.59499752 0.04774315 545.963385  42.4718352  43.1363939
LP       -30.6588754 0.79841066 0.01836798 518.072910 -33.0652295 -30.3900004
                  UB
beta[1]  -10.8925810
beta[2]    0.3269875
Deviance  48.4916661
LP       -30.0465224

Adaptive Metropolis

High thinning with 800, but nice mixing.

Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Iterations = 80000, Status = 2000, Thinning = 800, Algorithm = "AM",
    Specs = list(Adaptive = 200, Periodicity = 200))

Acceptance Rate: 0.26118
Algorithm: Adaptive Metropolis
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
     beta[1]      beta[2]
11.850286307  0.008006239

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
        All Stationary
Dbar 44.437     44.437
pD    1.630      1.630
DIC  46.067     46.067
Initial Values:
[1] -10   0

Iterations: 80000
Log(Marginal Likelihood): NA
Minutes of run-time: 0.36
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 0
Recommended Burn-In of Un-thinned Samples: 0
Recommended Thinning: 800
Specs: (NOT SHOWN HERE)
Status is displayed every 2000 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 100
Thinning: 800


Summary of All Samples
                Mean         SD        MCSE ESS          LB      Median
beta[1]  -11.7026349 1.97444818 0.188498237 100 -15.4468354 -11.6906840
beta[2]    0.2905268 0.04989136 0.004767928 100   0.2047257   0.2898906
Deviance  44.4368803 1.80543310 0.166862161 100  42.5326068  43.8582816
LP       -31.0345215 0.90945995 0.084411163 100 -33.4044687 -30.7391583
                  UB
beta[1]   -8.3034612
beta[2]    0.3801087
Deviance  49.1461781
LP       -30.0717992


Summary of Stationary Samples
                Mean         SD        MCSE ESS          LB      Median
beta[1]  -11.7026349 1.97444818 0.188498237 100 -15.4468354 -11.6906840
beta[2]    0.2905268 0.04989136 0.004767928 100   0.2047257   0.2898906
Deviance  44.4368803 1.80543310 0.166862161 100  42.5326068  43.8582816
LP       -31.0345215 0.90945995 0.084411163 100 -33.4044687 -30.7391583
                  UB
beta[1]   -8.3034612
beta[2]    0.3801087
Deviance  49.1461781
LP       -30.0717992

Adaptive Metropolis-within-Gibbs

A nice example where LaplacesDemon first makes me increase the number of samples then suddenly I get a run where it decides no burn in was needed. Needless to say with over a million samples I did not go to increase the thinning a bit more.

Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Iterations = 1200000, Status = 2000, Thinning = 41, Algorithm = "AMWG",
    Specs = list(Periodicity = 200))

Acceptance Rate: 0.24915
Algorithm: Adaptive Metropolis-within-Gibbs
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
 beta[1]  beta[2]
2.835066 2.835066

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
        All Stationary
Dbar 44.407     44.407
pD    1.896      1.896
DIC  46.304     46.304
Initial Values:
[1] -10   0

Iterations: 1200000
Log(Marginal Likelihood): NA
Minutes of run-time: 4.95
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 2
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 0
Recommended Burn-In of Un-thinned Samples: 0
Recommended Thinning: 44
Specs: (NOT SHOWN HERE)
Status is displayed every 2000 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 29268
Thinning: 41


Summary of All Samples
               Mean         SD        MCSE       ESS          LB     Median
beta[1]  -11.826585 1.96434369 0.146404922  39.11735 -15.8139146 -11.715120
beta[2]    0.293541 0.05072947 0.003808978  38.10170   0.1983134   0.290513
Deviance  44.407436 1.94752859 0.026815581 197.45010  42.4872692  43.822249
LP       -31.021258 0.98017230 0.013578417 192.69216 -33.7006842 -30.727256
                  UB
beta[1]   -8.0854510
beta[2]    0.3961088
Deviance  49.7428566
LP       -30.0525790


Summary of Stationary Samples
               Mean         SD        MCSE       ESS          LB     Median
beta[1]  -11.826585 1.96434369 0.146404922  39.11735 -15.8139146 -11.715120
beta[2]    0.293541 0.05072947 0.003808978  38.10170   0.1983134   0.290513
Deviance  44.407436 1.94752859 0.026815581 197.45010  42.4872692  43.822249
LP       -31.021258 0.98017230 0.013578417 192.69216 -33.7006842 -30.727256
                  UB
beta[1]   -8.0854510
beta[2]    0.3961088
Deviance  49.7428566
LP       -30.0525790

Common Code

All calls use the same data, model and plotting as given below.
s1 <- scan(what=list(integer(),double(),double()),text='
6 0 25.7  8 2 35.9 5 2 32.9 7 7 50.4 6 0 28.3 
7 2 32.3  5 1 33.2 8 3 40.9 6 0 36.5 6 1 36.5
6 6 49.6  6 3 39.8 6 4 43.6 6 1 34.1 7 1 37.4
8 2 35.2  6 6 51.3 5 3 42.5 7 0 31.3 3 2 40.6')
set2 <- data.frame(n=s1[[1]],x=s1[[3]],y=s1[[2]])
set2
glm(cbind(y,n-y)~ x,data=set2,family=binomial)
library(LaplacesDemon)
J <- 2 #Number of parameters
mon.names <- "LP"
parm.names <- c("beta[1]","beta[2]")
PGF <- function(Data) return(rnormv(Data$J,0,1000))
MyData <- list(J=J, PGF=PGF, n=set2$n, mon.names=mon.names,
    parm.names=parm.names, x=set2$x, y=set2$y)

Model <- function(parm, Data)
{
    ### Parameters
    beta <- parm[1:Data$J]
    ### Log-Prior
    beta.prior <- sum(dnormv(beta, 0, 1000, log=TRUE))
    ### Log-Likelihood
    mu <- beta[1] + beta[2]*Data$x
    p <- invlogit(mu)
    LL <- sum(dbinom(Data$y, Data$n, p, log=TRUE))
    ### Log-Posterior
    LP <- LL + beta.prior
    Modelout <- list(LP=LP, Dev=-2*LL, Monitor=LP,
        yhat=rbinom(length(p), Data$n, p), parm=parm)
    return(Modelout)
}
Initial.Values <- c(-10,0)


myplot <- function(where='none') {
    if (where!='none') {
        png(paste(where,'.png',sep=''))
    }
    par(mfrow=c(2,2),mar=rep(2,4),oma=c(0,0,2,0))
    plot(ts(Fit$Posterior1[,1]),ylab='')
    plot(ts(Fit$Posterior2[,1]),ylab='')
    plot(ts(Fit$Posterior1[,2]),ylab='')
    plot(ts(Fit$Posterior2[,2]),ylab='')
    title(main=Fit$Algorithm,outer=TRUE)
    if (where!='none') dev.off()
}

Sunday, October 5, 2014

Bayes models from SAS PROC MCMC in R, post 2

This is my second post in converting SAS's PROC MCMC examples in R. The task in his week is determining the transformation parameter in a Box-Cox transformation. SAS only determines Lambda, but I am not so sure about that. What I used to do was get an interval for Lambda, select an interpretable value (e.g. 1/2 is square root, -1 is inverse) and progress using this parameter. However, if that is my approach in a Bayesian context, maybe I should make my prior on lambda only those values I find acceptable. On the other hand, if I model lambda continuous, should not the uncertainty in Lambda also effect the uncertainty in the remainder of the model? In this post I will consider them all continuous. Hence the first step is to get SAS to give me those values, for which I added those parameters in the code and ran it through SAS Studio. Subsequently find the code for the various R approaches. For readability I stuck the code next to text, while output is at the bottom of the post. 

Note on previous post

Last week I failed to get some samplers. I now realize that some of these might have been solved with the ones trick or the zero trick. In addition, LaplacesDemon has dnormm() which is more elegant than my approach.

SAS

This is just the model from the SAS manual with extra parameters monitored.
proc mcmc data=boxcox nmc=50000 propcov=quanew seed=12567 monitor=(lda beta0 beta1 sd);
    ods select PostSumInt TADpanel;
    parms beta0 0 beta1 0 lda 1 s2 1;
    beginnodata;
    prior beta: ~ general(0);
    prior s2 ~ gamma(shape=3, scale=2);
    prior lda ~ unif(-2, 2);
    sd=sqrt(s2);
    endnodata;
    ys=(y**lda-1)/lda;
    mu=beta0+beta1*x;
    ll=(lda-1)*log(y)+lpdfnorm(ys, mu, sd);
    model general(ll);
run;

STAN

Rstan has the Jacobian needed for this exercise. It was fairly straightforward to fit it all together.
library(plyr)
library(dplyr)
# stan
library(rstan)
datain <- list(n=100,x=set2$x,y=set2$y)
parameters <- c('lda','beta0','beta1','std')
fit <- '
    data {
    int <lower=1> n;
    vector[n] y;
    vector[n] x;
    }
    parameters {
    real <lower=-2,upper=2> lda;
    real beta0;
    real beta1;
    real <lower=0> s2;
    }
    transformed parameters {
    vector[n] ys;
    vector[n] ll;
    vector[n] mu;
    real <lower=0> std;
    std <- sqrt(s2);
    for (i in 1:n) {
    ys[i] <- (pow(y[i],lda)-1)/lda;
    }
    mu <- beta0+beta1*x;
    ll <- (lda-1)*log(y);
    }
    model {
    s2 ~ gamma(3, 0.5);
    increment_log_prob(ll);
    increment_log_prob(normal_log(ys,mu,std));
    }
    ' %>%
stan(model_code = .,
    data = datain,
    pars=parameters)
fit

JAGS 

The JAGS code is not straightforward for two reasons.
The first reason is that Jags does not have the concept of Jacobian. However, upon some investigation I found the zeros trick. Which basically means that an arbitrary contribution to the deviance might be added by transforming that deviance and pretending to JAGS the values are from a Poisson distribution.
The second reason is that JAGS has directed acyclic graphs. At least, I think that is the reason. One cannot first calculate transformed values and then make them dependent on a distribution. My solution was to subtract transformed and estimated values and make the result normal around 0.
To get acceptable results in terms of n.eff and Rhat I increased the number of samples to 10000 per chain.
library(R2jags)
datain <- list(n=100,
    x=set2$x,
    y=set2$y,
    zero1 = rep(0,100),
    zero2 = rep(0,100))
parameters <- c('LDA','beta0','beta1','std')

jmodel <- function() {
    tau <- 1/s2
    std <- sqrt(s2)
    for (i in 1:n) {       
        mu[i] <- beta0+x[i]*beta1
        ys[i] <- (pow(y[i],LDA)-1)/LDA -mu[i]
        zero1[i] ~ dnorm(ys[i],tau)
        ll[i] <- -((LDA-1)*log(y[i]))+100
        zero2[i] ~ dpois(ll[i])
    }
    beta0 ~ dnorm(0,.0001)
    beta1 ~ dnorm(0,.0001)
    s2 ~ dgamma(3,.5)
    LDA ~ dunif(-2,2)
}
inits  <- function() {
    list(beta0=rnorm(1),
        beta1=rnorm(1),
        LDA=runif(1,-1,1),
        s2=runif(1,.1,3))
}
jj <- jags(model.file=jmodel,
    data=datain,
    inits=inits,
    n.iter=10000,
    parameters=parameters)
jj

MCMC

As such it is relatively easy to build the model for MCMC. Based on the demo vignette I was also able to set scale to obtain an acceptance rate about 0.2. After some additional calculations I came to batch length of 10000, and chose for 20 batches to keep computation time somewhat limited. All in all it took quite some attempts to get to these settings in terms of nbatch and batch length. MCMC only gives summary statistics based on each separate batch. These are subsequently summarized. The results seem a bit different from SAS and JAGS. In addition a small tweak was needed to get SD as output rather than internal variable S2.
mydens <- function(x,...) {
    if (x[4]<0 ) return(-Inf)
    std <- sqrt(x[4])
    mu <- x[2]+set2$x*x[3]
    ys <- (set2$y^x[1]-1)/x[1]
    ll <- (x[1]-1)*log(set2$y)
    sum(dnorm(ys,mu,std,log=TRUE)+ll) +dgamma(x[4],shape=3,scale=2)
}
sam <-  metrop(mydens,
    initial=c(.46,-2,1.8,1.3),
    nbatch=20,
    blen=10000,
    scale=.025,
    outfun=function(z, ...) {
        z[4] <- sqrt(z[4])
        c(z, z^2)}
)
sam$accept
foo <- apply(sam$batch, 2, mean)
mu <- foo[1:4]
sigmasq <- foo[5:8] - mu^2
data.frame(param=c('LDA','beta0','beta1','SD'),mu,sigmasq)



MCMCpack

Obviously, the deviance code is same as MCMC. Based on the ACF a thinning of 20 was chosen. I did not see an obvious reason to change from the default 20000 samples.
library(MCMCpack)
library(coda)
mydens <- function(x,...) {
    if (x[4]<0 ) return(-Inf)
    std <- sqrt(x[4])
    mu <- x[2]+set2$x*x[3]
    ys <- (set2$y^x[1]-1)/x[1]
    ll <- (x[1]-1)*log(set2$y)
    sum(dnorm(ys,mu,std,log=TRUE)+ll) +dgamma(x[4],shape=3,scale=2)
}
sam <- MCMCmetrop1R(mydens,
    theta.init=c(.46,-2,1.8,1.3),
    tune=1.5,thin=20)
summary(sam)

rcppbugs

Rcppbugs has many of the same limits as JAGS in terms of available distributions. The same tricks as used in JAGS are hence used in rcppbugs, with one exception, where JAGS has a Poison distribution, rcppbugs lacks it, hence the ones trick is used, which is the equivalent to zeros trick, but using Bernoulli distribution. The ACF gave quite a big lag, hence I upped both the number of samples and increased thinning.

library(rcppbugs)
x <- set2$x
y <- set2$y
one <- rep(1,100)
yzero <- rep(0,100)
beta <- mcmc.normal(c(.1,.1),mu=0,tau=0.0001)
LDA <- mcmc.uniform(.5,-2,2)
tau.y <- mcmc.gamma(sd(as.vector(y)),alpha=0.1,beta=0.1)
y.hat <- deterministic(function(x,beta,y,LDA) {
        ypred=beta[1]+beta[2]*x
        ytrans=(y^LDA-1)/LDA
        ypred-ytrans
    }, x, beta,y,LDA)
Jac <- deterministic(function(LDA,y)
        exp((LDA-1)*log(y))/100,LDA,y)
cJac <- mcmc.bernoulli(one,Jac,observed=TRUE)
y.lik <- mcmc.normal(yzero,mu=y.hat,tau=tau.y,observed=TRUE)

m <- create.model(beta,LDA, tau.y,
    y.hat,
    y.lik,
    Jac,
    cJac)
ans <- run.model(m, iterations=1e6L, burn=1e4L, adapt=1e3L, thin=1000L)

mysum <- function(x) c(mean=mean(x),sd=sd(x),quantile(x,c(.025,.5,.975)))
rbind(LDA=mysum(ans$LDA),
    beta0=mysum(ans$b[,1]),
    beta1=mysum(ans$b[,2]),
    sd=mysum(1/sqrt(ans$tau.y)))

LaplacesDemon

LaplacesDemon also needed a lot of samples. At a later attempt I tried the HARM sampler which worked much better. Getting the right sampler is important.
library('LaplacesDemon')
mon.names <- "LP"
parm.names <- c('beta0','beta1','LDA','std')
PGF <- function(Data) c(rnorm(2,0,1),runif(1,.4,.5),runif(1,.5,2))

MyData <- list(mon.names=mon.names,
    parm.names=parm.names,
    PGF=PGF,
    y=set2$y,
    x=set2$x)
Model <- function(parm, Data)
{
    ### Parameters
    std <- interval(parm[4], 1e-100, Inf)
    parm[4] <- std
    s2 <- parm[4]^2
    mu <- parm[1]+Data$x*parm[2]
    ys <- (Data$y^parm[3]-1)/parm[3]
    ll <- (parm[3]-1)*log(Data$y)
    LL <- sum(dnorm(ys,mu,std,log=TRUE)+ll)
    LP <- LL+dgamma(s2,shape=3,scale=2)
    Modelout <- list(LP=LP, Dev=-2*LL, Monitor=LP,
        yhat=1,
        parm=parm)
    return(Modelout)
}
# very dependent on initial values
Initial.Values <- GIV(Model, MyData, PGF=TRUE)
Fit <- LaplacesDemon(Model,
    Data=MyData,
    Algorithm='HARM',
    Specs=list(alpha.star=0.234, B=NULL),
    Initial.Values = Initial.Values,
    Iterations=40000,Status=2000,Thinning=30
)
Fit

Output

SAS


                              Standard 
Parameter     N     Mean      Deviation     95% HPD Interval
lda        50000    0.4701     0.0290     0.4142     0.5281
beta0      50000   -2.2064     0.4215    -3.0564    -1.4201
beta1      50000    1.9260     0.1806     1.5854     2.2826
sd         50000    1.3275     0.1482     1.0648     1.6328


STAN

Inference for Stan model: __LHS.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.

         mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
lda      0.47    0.00 0.03    0.41    0.45    0.47    0.49    0.53   595 1.01
beta0   -2.20    0.02 0.40   -3.01   -2.45   -2.18   -1.94   -1.44   684 1.01
beta1    1.92    0.01 0.17    1.59    1.81    1.91    2.02    2.28   571 1.01
std      1.32    0.01 0.15    1.06    1.22    1.31    1.42    1.65   691 1.00
lp__  -309.43    0.05 1.45 -313.02 -310.13 -309.09 -308.36 -307.64   977 1.00

Samples were drawn using NUTS(diag_e) at Sat Oct  4 15:14:25 2014.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).

JAGS 


Inference for Bugs model at "/tmp/RtmpSgjWNR/modelef330519173.txt", fit using jags,
 3 chains, each with 10000 iterations (first 5000 discarded), n.thin = 5
 n.sims = 3000 iterations saved
           mu.vect sd.vect      2.5%       25%       50%       75%     97.5%
LDA          0.472   0.030     0.413     0.451     0.472     0.491     0.535
beta0       -2.229   0.443    -3.200    -2.512    -2.199    -1.920    -1.450
beta1        1.938   0.192     1.608     1.805     1.925     2.054     2.356
std          1.336   0.155     1.079     1.227     1.319     1.433     1.680
deviance 20620.595   3.365 20616.435 20618.142 20619.763 20622.197 20629.012
          Rhat n.eff
LDA      1.009  2900
beta0    1.003  1400
beta1    1.008  1500
std      1.002  3000
deviance 1.011   480

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)
pD = 5.6 and DIC = 20626.2
DIC is an estimate of expected predictive error (lower deviance is better).

MCMC


  param         mu      sigmasq
1   LDA  0.4614172 0.0007537545
2 beta0 -2.0808097 0.1414263712
3 beta1  1.8693639 0.0267559710
4    SD  1.2832247 0.0180115043

MCMCpack

Note that MCMCpack has s2 rather than sd in output

Iterations = 501:20481
Thinning interval = 20
Number of chains = 1
Sample size per chain = 1000

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

        Mean      SD  Naive SE Time-series SE
[1,]  0.4657 0.02927 0.0009256      0.0009256
[2,] -2.1565 0.40274 0.0127358      0.0113576
[3,]  1.8988 0.17764 0.0056175      0.0052342
[4,]  1.7048 0.36930 0.0116782      0.0133882

2. Quantiles for each variable:

        2.5%     25%    50%     75%   97.5%
var1  0.4058  0.4479  0.466  0.4854  0.5215
var2 -2.9944 -2.4064 -2.146 -1.8799 -1.4346
var3  1.5606  1.7810  1.890  2.0099  2.2647
var4  1.1204  1.4473  1.652  1.9261  2.5331

rcppbugs


            mean         sd       2.5%        50%      97.5%
LDA    0.4597165 0.02875052  0.4050122  0.4592264  0.5207108
beta0 -2.0937414 0.39740763 -2.9511520 -2.0668905 -1.4015904
beta1  1.8647679 0.17411536  1.5640795  1.8566258  2.2509009
sd     1.2612126 0.14011464  1.0228554  1.2502009  1.5634482

LaplacesDemon

Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Iterations = 40000, Status = 2000, Thinning = 30, Algorithm = "HARM",
    Specs = list(alpha.star = 0.234, B = NULL))

Acceptance Rate: 0.22825
Algorithm: Hit-And-Run Metropolis
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
       beta0        beta1          LDA          std
0.1329756466 0.0283948376 0.0007961449 0.0174295808

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
         All Stationary
Dbar 619.952    619.768
pD     9.531      3.530
DIC  629.483    623.298
Initial Values:
[1] -0.2820944  0.3190142  0.4778030  1.3380326

Iterations: 40000
Log(Marginal Likelihood): -313.4286
Minutes of run-time: 0.12
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 4
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 532
Recommended Burn-In of Un-thinned Samples: 15960
Recommended Thinning: 31
Specs: (NOT SHOWN HERE)
Status is displayed every 2000 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 1333
Thinning: 30


Summary of All Samples
                 Mean        SD        MCSE       ESS           LB       Median
beta0      -2.0923727 0.3614098 0.091343057  17.08625   -2.8552455   -2.0861434
beta1       1.8711745 0.1631220 0.037952180  22.13164    1.5969508    1.8620763
LDA         0.4612769 0.0282230 0.004279672  34.00465    0.4086439    0.4616171
std         1.2653917 0.1320557 0.020805981  98.82246    1.0352842    1.2574597
Deviance  619.9523160 4.3660058 0.415674235 188.93047  616.3618665  619.2464350
LP       -309.9042345 2.1778219 0.206890192 190.00073 -313.2903260 -309.5529196
                   UB
beta0      -1.4732055
beta1       2.2047884
LDA         0.5171681
std         1.5495557
Deviance  626.7398118
LP       -308.1117717


Summary of Stationary Samples
                 Mean         SD        MCSE        ESS           LB
beta0      -2.0498690 0.36586611 0.118090993   8.836406   -2.9348999
beta1       1.8509120 0.16921671 0.050214255  15.798638    1.5936363
LDA         0.4577846 0.02886563 0.006775378  21.788913    0.4078682
std         1.2449204 0.12994746 0.027000671  55.838883    1.0332563
Deviance  619.7681778 2.65715568 0.326080727 121.032423  616.3311632
LP       -309.8148748 1.32309998 0.206890192 123.122008 -313.0799005
               Median           UB
beta0      -1.9902222   -1.5178152
beta1       1.8236876    2.2326418
LDA         0.4558953    0.5206561
std         1.2320277    1.5344502
Deviance  619.3065650  626.3929182
LP       -309.5806926 -308.1020531

Appendix: data used


s1 <- scan(what=list(double(),double()),text='
10.0  3.0  72.6  8.3  59.7  8.1  20.1  4.8  90.1  9.8   1.1  0.9
78.2  8.5  87.4  9.0   9.5  3.4   0.1  1.4   0.1  1.1  42.5  5.1
57.0  7.5   9.9  1.9   0.5  1.0 121.1  9.9  37.5  5.9  49.5  6.7
8.3  1.8   0.6  1.8  53.0  6.7 112.8 10.0  40.7  6.4   5.1  2.4
73.3  9.5 122.4  9.9  87.2  9.4 121.2  9.9  23.1  4.3   7.1  3.5
12.4  3.3   5.6  2.7 113.0  9.6 110.5 10.0   3.1  1.5  52.4  7.9
80.4  8.1   0.6  1.6 115.1  9.1  15.9  3.1  56.5  7.3  85.4  9.8
32.5  5.8  43.0  6.2   0.1  0.8  21.8  5.2  15.2  3.5   5.2  3.0
0.2  0.8  73.5  8.2   4.9  3.2   0.2  0.3  69.0  9.2   3.6  3.5
0.2  0.9 101.3  9.9  10.0  3.7  16.9  3.0  11.2  5.0   0.2  0.4
80.8  9.4  24.9  5.7 113.5  9.7   6.2  2.1  12.5  3.2   4.8  1.8
80.1  8.3  26.4  4.8  13.4  3.8  99.8  9.7  44.1  6.2  15.3  3.8
2.2  1.5  10.3  2.7  13.8  4.7  38.6  4.5  79.1  9.8  33.6  5.8
9.1  4.5  89.3  9.1   5.5  2.6  20.0  4.8   2.9  2.9  82.9  8.4
7.0  3.5  14.5  2.9  16.0  3.7  29.3  6.1  48.9  6.3   1.6  1.9
34.7  6.2  33.5  6.5  26.0  5.6  12.7  3.1   0.1  0.3  15.4  4.2
2.6  1.8  58.6  7.9  81.2  8.1  37.2  6.9
')
set2 <- data.frame(x=s1[[2]],y=s1[[1]])