// 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); // calculate the posterior parameters a_tilde = N + 1; b_tilde = sumy + 0.25; // draw samples from the posterior G = 100; x = gamrnd(a_tilde, b_tilde, G, 1); // calculate the moments and the Monte Carlo standard error E_lambda = mean(x); V_lambda = var(x); MCse = sqrt(V_lambda/G); // print the results print( [E_lambda, V_lambda, MCse] );