// import the data and get N and the sum of the values in y Data = webimport("www.bayeconsoft.com/datasets/WaitingTimes.csv"); N = rows(Data); sumy = sum(Data.y); // set the values of the hyperparameters and the tuning parameter m = 1.4; s2 = 0.9; tuning = 0.3; // set the number of iterations D = 3000; // # of burn-in iterations G = 10000; // # of retained draws // set the starting value for lambda and calculate its logarithm lambda = 2; loglambda = log(lambda); // initialize a vector to store the draws draws = zeros(G,1); // start the algorithm for (g=1:D+G) // draw log-lambda star from the proposal and calculate its exponential loglambda_star = normrnd(loglambda,tuning); lambda_star = exp(loglambda_star); // calculate the logarithm of the Metropolis-Hastings ratio logMH = N*(loglambda_star-loglambda) - (lambda_star-lambda)*sumy - ((loglambda_star-m)^2 - (loglambda-m)^2)/(2*s2); // accept/reject the proposed move if ( log(unifrnd()) < logMH ) lambda = lambda_star; loglambda = loglambda_star; end // store the results from the current iteration ================ if (g>D) draws(g-D) = lambda; end end // summarize the draws from the posterior print( [mean(draws); var(draws)] );