{"word_count": 426, "character_count": 3369}
This afternoon, while I was discussing with Montserrat (aka @mguillen_estany) we were wondering how long it might take to run a regression model. More specifically, how long it might take if we use a Bayesian approach. My guess was that the time should probably be linear in , the number of observations. But I thought I would be good to check.Let us generate a big dataset, with one million rows,&#62; n=1e6&#62; X=runif(n)&#62; Y=2+5*X+rnorm(n)&#62; B=data.frame(X,Y)Consider as a benchmark the standard linear regression,&#62; lm_freq = function(n){+   idx = sample(1:1e6,size=n)+   reg = lm(Y~X,data=B[idx,])+   summary(reg)+ }Here the regression is a subset of smaller size. We can do the same with a Bayesian approach, using stan,&#62; stan_lm ="+ data {+ int N;+ vector[N] x;+ vector[N] y;+ }+ parameters {+ real alpha;+ real beta;+ real tau;+ }+ transformed parameters {+ real sigma;+ sigma &#60;- 1 / sqrt(tau);+ }+ model{+ y ~ normal(alpha + beta * x, sigma);+ alpha ~ normal(0, 10);+ beta ~ normal(0, 10);+ tau ~ gamma(0.001, 0.001);+ }+ "Define then the model&#62; library(rstan)&#62; system.time(   stanmodel &#60;&#60;- stan_model(model_code = stan_lm))utilisateur     système      écoulé       0.043       0.000       0.043We want to see how long it might take to run a regression,&#62; lm_bayes = function(n){+   idx = sample(1:1e6,size=n)+   fit = sampling(stanmodel,+       data = list(N=n,+                   x=X[idx],+                   y=Y[idx]),+       iter = 1000, warmup=200)+   summary(fit)+ }We use the following package to see how long it takes&#62; library(microbenchmark)&#62; time_lm = function(n){+  M = microbenchmark(lm_freq(n),+      lm_bayes(n),times=50)+  return(apply( matrix(M$time,nrow=2),1,mean))+ }We can now compare the time it took with ten, one hundred, on thousand, and ten thousand observations,&#62; vN = c(10,100,1000,10000)&#62; T = Vectorize(time_lm)(vN)we can then plot it&#62; plot(vN,T[2,]/1e6,log="xy",col="red",type="b",+      xlab="Number of Observations",ylab="Time")&#62; lines(vN,T[1,]/1e6,col="blue",type="b")It looks like (if we forget about the very small sample) that the time it takes to run a regression is linear, with the two techniques (the frequentist and the Bayesian ones).And actually, the same story olds for logistic regressions. Consider the following dataset&#62; n=1e6&#62; X=runif(n)&#62; S=-3+2*X+rnorm(n)&#62; Y=rbinom(n,size=1,prob=exp(S)/(1+exp(S)))&#62; B=data.frame(X,Y)The frequentist version of the logistic regression is&#62; glm_freq = function(n){+   idx = sample(1:1e6,size=n)+   reg = glm(Y~X,data=B[idx,],family=binomial)+   summary(reg)+ }and the Bayesian one, using stan,&#62; stan_glm = "+ data {+ int N;+ vector[N] x;+ int&#60;lower=0,upper=1&#62; y[N];+ }+ parameters {+ real alpha;+ real beta;+ }+ model {+ alpha ~ normal(0, 10);+ beta ~ normal(0, 10);+ y ~ bernoulli_logit(alpha + beta * x);+ }+ "&#62; stanmodel = stan_model(model_code = stan_glm) )&#62; glm_bayes = function(n){+   idx = sample(1:1e6,size=n)+   fit = sampling(stanmodel,+        data = list(N=n,+        x = X[idx],+        y = Y[idx]),+        iter = 1000, warmup=200)+   summary(fit)+ }Again, we can see how long it takes to run those regression models&#62; time_gl m= function(n){+   M = microbenchmark(glm_freq(n),+   glm_bayes(n),times=50)+   return(apply( matrix(M$time,nrow=2),1,mean))+ }&#160;