Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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 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  beta
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:
 -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   -9.0326306 1.04621897 0.220388172  7.322802 -11.2499597  -8.9042795
beta    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   -7.1874098
beta    0.2772263
Deviance  50.1501559
LP       -30.0892713

Summary of Stationary Samples
Mean         SD        MCSE      ESS          LB      Median
beta   -8.9359418 0.99836162 0.349118267 3.269911 -10.6298830  -8.8780160
beta    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   -7.4135378
beta    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     beta
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:
 -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  -12.5516619 1.69184474 0.251682993  43.36212 -15.7525684 -12.537176
beta    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   -9.3734152
beta    0.3968852
Deviance  49.3254872
LP       -30.0617562

Summary of Stationary Samples
Mean         SD        MCSE       ESS          LB     Median
beta  -12.5516619 1.69184474 0.251682993  43.36212 -15.7525684 -12.537176
beta    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   -9.3734152
beta    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     beta
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:
 -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  -11.8540232 1.90772737 0.141469541  28.83441 -15.618594 -11.8276283
beta    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   -8.489146
beta    0.392589
Deviance  49.076360
LP       -30.058998

Summary of Stationary Samples
Mean         SD        MCSE       ESS          LB      Median
beta  -11.7694571 1.87121305 0.153383946  24.00847 -15.4323928 -11.8230837
beta    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   -8.3411645
beta    0.3869563
Deviance  48.9102683
LP       -30.0612554

#### 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
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
beta     beta
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:
 -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  -11.5514694 1.92227432 0.252091958  40.70347 -15.7178588 -11.4718511
beta    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   -8.0953413
beta    0.3936748
Deviance  49.2894939
LP       -30.0514927

Summary of Stationary Samples
Mean         SD        MCSE       ESS          LB      Median
beta  -11.5514694 1.92227432 0.252091958  40.70347 -15.7178588 -11.4718511
beta    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   -8.0953413
beta    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. 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
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
beta     beta
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:
 -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  -11.314847 2.17303791 0.51002442 11.27064 -15.1041582 -11.4860852
beta    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   -7.2420670
beta    0.3795564
Deviance  49.4135552
LP       -30.0677263

Summary of Stationary Samples
Mean         SD        MCSE       ESS          LB      Median
beta  -13.3566117 0.72016481 0.299751317  5.473664 -15.0713103 -13.1769958
beta    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  -12.2282779
beta    0.3766297
Deviance  47.5762853
LP       -30.2209601

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 beta 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:  -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 -10.922145 4.4338567 0.50560810 80 -16.5329510 -11.2359881 beta 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 -5.6863543 beta 0.4200815 Deviance 62.5858557 LP -30.0729465 Summary of Stationary Samples Mean SD MCSE ESS LB Median beta -10.922145 4.4338567 0.50560810 80 -16.5329510 -11.2359881 beta 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 -5.6863543 beta 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 beta 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:  -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 -11.4418192 1.06961398 0.191704524 7.357239 -14.1113189 -11.4442594 beta 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 -9.6350028 beta 0.3515657 Deviance 47.9431739 LP -30.0443227 Summary of Stationary Samples Mean SD MCSE ESS LB Median beta -11.9854873 0.56317415 0.15499286 8.169076 -12.9867511 -11.9369396 beta 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 -10.8925810 beta 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 beta 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:  -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 -11.7026349 1.97444818 0.188498237 100 -15.4468354 -11.6906840 beta 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 -8.3034612 beta 0.3801087 Deviance 49.1461781 LP -30.0717992 Summary of Stationary Samples Mean SD MCSE ESS LB Median beta -11.7026349 1.97444818 0.188498237 100 -15.4468354 -11.6906840 beta 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 -8.3034612 beta 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 beta 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:  -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 -11.826585 1.96434369 0.146404922 39.11735 -15.8139146 -11.715120 beta 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 -8.0854510 beta 0.3961088 Deviance 49.7428566 LP -30.0525790 Summary of Stationary Samples Mean SD MCSE ESS LB Median beta -11.826585 1.96434369 0.146404922 39.11735 -15.8139146 -11.715120 beta 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 -8.0854510 beta 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[],x=s1[],y=s1[]) 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","beta") 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 + beta*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()
}