Speed comparison between R, Julia and Fortran on Bayesian probit regression (see: https://taptoenix.wordpress.com/2015/10/31/speed-comparison-between-r-julia-and-fortran-on-bayesian-probit-regression/)

readme

R, Julia and Fortran code used for speed comparison on Bayesian probit regression. 
See blogpost for more details: https://taptoenix.wordpress.com/2015/10/31/speed-comparison-between-r-julia-and-fortran-on-bayesian-probit-regression/

binprobbayes.f90

! Written by Dries F. Benoit 
! Faculty of economics and business administration
! Ghent University - BELGIUM

! Draws from the posterior distribution of the Bayesian probit model
! with normal prior on the regression coefficients.

! Algorithm based on:
! Albert and Chib (1993). Bayesian Analysis of Binary and Polychotomous
! Response Data. Journal of the American Statistical Association, 
! 88(422), 669-679.

! Licence GLPv3 (or higher)

subroutine binprobbayes(y,x,b0,bb0,r,thin,n,k,betadraw)

implicit none

! Precision statement
integer, parameter :: dp = kind(1.0d0)

! Input arguments
integer, intent(in) :: r,thin,n,k
integer, dimension(n), intent(in) :: y 
real(dp), dimension(n,k), intent(in) :: x 
real(dp), dimension(k), intent(in) :: b0 
real(dp), dimension(k,k), intent(in) :: bb0 

! Output arguments
real(dp), dimension(r,k), intent(out) :: betadraw 

! Internal arguments
integer :: ok,i,ii
real(dp), dimension(k) :: b,beta
real(dp), dimension(k,k) :: bb
real(dp), dimension(n) :: ystar 

! Set starting values
beta = 0.0_dp
ystar = 0.0_dp

bb = bb0
call dpotrf('U',k,bb,k,ok)
call dpotri('U',k,bb,k,ok)
do i = 1,(k-1)
        do ii = (i+1),k
                bb(ii,i) = bb(i,ii)
        enddo
enddo

b = matmul(bb,b0)
bb = bb + matmul(transpose(x),x)

call dpotrf('U',k,bb,k,ok)
call dpotri('U',k,bb,k,ok)
do i = 1,(k-1)
        do ii = (i+1),k
                bb(ii,i) = bb(i,ii)
        enddo
enddo

! Start mcmc
do i = 1,r
        
        ! Draw new value for ystar
        do ii = 1,n
                if (y(ii) .eq. 0) then
                        call rtnorm(0.0_dp,.FALSE.,dot_product(x(ii,:),beta),1.0_dp,ystar(ii))
                else
                        call rtnorm(0.0_dp,.TRUE.,dot_product(x(ii,:),beta),1.0_dp,ystar(ii))
                endif
        enddo

        ! Draw new value for beta
        call rmvnorm(matmul(bb,(b+matmul(transpose(x),ystar))),bb,k,beta)

        ! Save draw
        betadraw(i,:) = beta

enddo

end subroutine binprobbayes 

binprobbayes.jl

# Written by Dries F. Benoit 
# Faculty of economics and business administration
# Ghent University - BELGIUM

# Draws from the posterior distribution of the Bayesian probit model
# with normal prior on the regression coefficients.

# Algorithm based on:
# Albert and Chib (1993). Bayesian Analysis of Binary and Polychotomous
# Response Data. Journal of the American Statistical Association, 
# 88(422), 669-679.

# Licence GLPv3 (or higher)

using Distributions
function binprobbayes(y::Array{Int64,1}, X::Array{Float64,2}, # Data
			b0::Array{Float64,1}, B0::Array{Float64,2}, # Prior
			R::Int64, thin::Int64) # Mcmc

  n = size(X,1)
  k = size(X,2)

  #B = invvarcov(invvarcov(B0) + *(X',X))
  #Bb = *(invvarcov(B0),b0)
  B = inv(inv(B0) + *(X',X))
  Bb = *(inv(B0),b0)

  # Initialize space to save draws
  #betadraw = reshape(repeat([NaN],inner=[R*k]),R,k)
  betadraw = Array(Float64,R,k)

  # Set starting values
  beta = repeat([0.],inner=[k])
  ystar = repeat([0.],inner=[n])
  
  # Start mcmc
  for i = 1:R

        # Draw new value for ystar
        for ii = 1:n
                if y[ii]==0
                        ystar[ii] = rtnorm(*(X[ii,:],beta)[1],1.,max=0.)
                else
                        ystar[ii] = rtnorm(*(X[ii,:],beta)[1],1.,min=0.)
                end
        end

        # Draw new value for beta
        beta = rand(MvNormal(B*(Bb+*(X',ystar)),B))

	# Save draw
        betadraw[i,:] = beta

  end

  return betadraw

end

binprobbayes.r

# Written by Dries F. Benoit 
# Faculty of economics and business administration
# Ghent University - BELGIUM

# Draws from the posterior distribution of the Bayesian probit model
# with normal prior on the regression coefficients.

# Algorithm based on:
# Albert and Chib (1993). Bayesian Analysis of Binary and Polychotomous
# Response Data. Journal of the American Statistical Association, 
# 88(422), 669-679.

# Licence GLPv3 (or higher)

binprobbayes <- function(y,X,b0,B0,R,thin){

	n <- nrow(X) 
	k <- ncol(X)

	B <- solve(solve(B0) + t(X)%*%X)
	Bb <- solve(B0)%*%b0

	# Initialize space to save draws
	betadraw <- matrix(NA,nrow=R,ncol=k)

	# Set starting values
	beta <- rep(0, k)
	ystar <- rep(0,n)

	# Start mcmc
	for (i in 1:R){
	
		# Draw next value for ystar
		for (ii in 1:n){
			if (y[ii]==0){
				ystar[ii] <- rtnorm(X[ii,]%*%beta,1,max=0)
			} else {
				ystar[ii] <- rtnorm(X[ii,]%*%beta,1,min=0)
			}
		}

		# Draw new value for beta
		beta <- B%*%(Bb + (t(X)%*%ystar)) + t(chol(B))%*%rnorm(n=k)

		# Save draw
		betadraw[i,] = beta
	}

	return(betadraw)
}

datagen.jl

# Generate data

## Number of observations 
n = 5000;

## Parameters
beta = [-1.; 1.; .5];
sigsq = 1.;

## Number of regression coefs
k = size(beta,1);

## Generate independent vars
X = rand(n,k-1);
X = [ones(n) X]; # concat column of ones

## Generate dependent var
y_star = X*beta + rand(Normal(0.,sigsq),n);
y = map((x) -> ifelse(x<0.,0,1),y_star);


# Save as csv

## First create dataframe
using DataArrays, DataFrames
df = hcat(DataFrame(y=y),convert(DataFrame,X));

## Write to disk
writetable("probitdata.csv", df)

probitdata.csv

rmvnorm.f90

! Written by Dries F. Benoit 
! Faculty of economics and business administration
! Ghent University - BELGIUM

! This code generates one draw from the multivariate normal 
! distribution. 
! This code is based on the Cholesky decomposition of the
! variance-covariance matrix.

! Input arguments:
!	- mean   : mean vector
!	- varcov : varcov matrix 
!	- k      : dimension of mvnorm 

! Output arguments:
!	- fn_val	: random draw from N(0,1) distribution

! Licence: GPLv3 (or higher) 

subroutine rmvnorm(mean,varcov,k,fn_val)

implicit none

! Precision statement
integer, parameter :: dp = kind(1.0d0)

! Input arguments
integer, intent(in) :: k
real(dp), dimension(k), intent(in) :: mean
real(dp), dimension(k,k), intent(in) :: varcov

! Output arguments
real(dp), dimension(k), intent(out) :: fn_val 

! Internal arguments
integer :: ok,i1,i2
real(dp), dimension(k,k) :: cvarcov


! Cholesky root of varcov
cvarcov = varcov
call dpotrf('U',k,cvarcov,k,ok)
do i1 = 1,(k-1)
  do i2 = (i1+1),k
    cvarcov(i2,i1) = 0.0_dp
  enddo
enddo

! Standard normal vector
do i1 = 1,k
  call rnorm(fn_val(i1))
enddo

! Linear transformation
fn_val = mean + matmul(fn_val,cvarcov)

end subroutine rmvnorm

rnorm.f90

! Written by Dries F. Benoit 
! Faculty of economics and business administration
! Ghent University - BELGIUM

! This code generates one draw from the standard normal 
! distribution. Note that more efficient code is possible
! when more than one normal draw is required.
! This code is based on the Box-Muller method.

! Output arguments:
!	- fn_val	: random draw from N(0,1) distribution

! Licence: GPLv3 (or higher) 


subroutine rnorm(fn_val)

implicit none

! precision statement:
integer, parameter :: dp = kind(1.0d0)

! output arguments:
real(dp), intent(out) :: fn_val

! internal arguments:
real(dp) :: pi
real(dp), dimension(1:2) :: u

pi = 3.14159265358979323846_dp

call random_number(u)

fn_val = sqrt(-2.0_dp*log(u(1))) * cos(2.0_dp*pi*u(2))

end subroutine rnorm

rtnorm.f90

! Written by Dries F. Benoit 
! Faculty of economics and business administration
! Ghent University - BELGIUM

! Returns one draw from the truncated normal distribution

! Algorithm based on:
! Geweke, J. (1991). Efficient Simulation From the Multivariate Normal 
! and Student t-Distributions Subject to Linear Constraints, in Computer 
! Sciences and Statistics Proceedings of the 23d Symposium on the 
! Interface, pp. 571-578.

! This subroutine makes use of the subroutines:
!	- rnorm		: Box-Muller method for random normal draws

! Input arguments:
! a             -	trucation point
! lb	        -	logical:        if .TRUE. then trucation (a,+Inf)
!		        		if .FALSE. then truncation (-Inf,a)
! mu	        -	mean of trunc normal
! sigma         -	sd of trunc normal
! fn_val        -	random draw from trunc normal

! Licence: GPLv3 (or higher)


subroutine rtnorm(a, lb, mu, sigma, fn_val)
  
implicit none

! Precision statement:
integer, parameter :: dp = kind(1.0d0)

! Input arguments:
logical, intent(in) :: lb
real(dp), intent(in) :: a, mu, sigma

! Output arguments:
real(dp), intent(out) :: fn_val

! Internal arguments:
real(dp) :: z, phi_z, az, c
real(dp), dimension(1:2) :: u


! Rescale truncation point
az=(a-mu)/sigma

if (lb) then
    c=az
  else 
    c=-az
endif

if (c<.45_dp) then

  ! normal rejection sampling
  do
    call rnorm(u(1))
    
    if (u(1)>c) exit
  end do
  z=u(1)

else

  ! exponential rejection sampling
  do
    ! Create exponential random variate z
    ! from uniform random variate u(1)
    call random_number(u)
    z = -log(u(1))/c

    phi_z = exp(-.5_dp * z**2_dp) !see Geweke
    if (u(2)<phi_z) exit
  end do
  z=z+c

end if

if (lb) then 
  fn_val = mu + sigma*z
else
  fn_val = mu - sigma*z
end if

end subroutine rtnorm

rtnorm.jl

# Written by Dries F. Benoit 
# Faculty of economics and business administration
# Ghent University - BELGIUM

# Returns one draw from the truncated normal distribution

# Algorithm based on:
# Geweke, J. (1991). Efficient Simulation From the Multivariate Normal 
# and Student t-Distributions Subject to Linear Constraints, in Computer 
# Sciences and Statistics Proceedings of the 23d Symposium on the 
# Interface, pp. 571-578.

# Input arguments:
# mu	        -	mean of trunc normal
# sigma         -	sd of trunc normal
# min           -       truncation point from below
# max           -       truncation point from above
# NOTE: truncation from both above and below is not supported

# Licence: GPLv3 (or higher)

using Distributions
function rtnorm(mu::Float64, sd::Float64; min::Float64=-Inf, max::Float64=Inf)

  local x

  if isfinite(max)
    lt = false
    c = -(max-mu)/sd
  else
    lt = true
    c = (min-mu)/sd
  end

  if c <= .45 
    # normal rejection sampling
    while true
      x = rand(Normal(0.,1.))
      if x > c
        break
      end
    end
  
  else
    # exponential rejection sampling
    while true
      x = rand(Exponential(1./c)) # check parametrization!
      u = rand()
      if u < exp(-.5*(x^2))
        x = x + c
        break
      end
    end
  end

  if lt
    return mu + sd * x
  else
    return mu - sd * x
  end
  
end

rtnorm.r

# Written by Dries F. Benoit 
# Faculty of economics and business administration
# Ghent University - BELGIUM

# Returns one draw from the truncated normal distribution

# Algorithm based on:
# Geweke, J. (1991). Efficient Simulation From the Multivariate Normal 
# and Student t-Distributions Subject to Linear Constraints, in Computer 
# Sciences and Statistics Proceedings of the 23d Symposium on the 
# Interface, pp. 571-578.

# Input arguments:
# mu	        -	mean of trunc normal
# sigma         -	sd of trunc normal
# min           -       truncation point from below
# max           -       truncation point from above
# NOTE: truncation from both above and below is not supported

# Licence: GPLv3 (or higher)

rtnorm <- function(mu, sd, min=-Inf, max=Inf){

  if(is.finite(max)){
    lt <- FALSE
    c <- -(max-mu)/sd
  }else{
    lt=TRUE
    c <- (min-mu)/sd
  }	
  
  if(c <= .45){
    # normal rejection sampling
    repeat{
      x <- rnorm(n=1, mean=0, sd=1)
      if(x>c)break
    }
  } else if (c > .45){
    # exponential rejection sampling
    repeat{
      x <- rexp(n=1, rate=c)
      u <- runif(n=1, min=0, max=1)
      if(u < exp(-.5*(x^2)))break
    }
    x <- x + c
  }

  if(lt){
    return(mu+sd*x)
  }else{
    return(mu-sd*x)
  }
}

test_bayesprob.r

# Source the algorithms
source("binprobbayes.r")
source("rtnorm.r")

# Read in the data
df <- read.csv("probitdata.csv")

# Execute algorithm
tic <- proc.time()
out <- binprobbayes(df$y,as.matrix(df[,2:4]),c(0,0,0),diag(3)*1000,5000,1)
proc.time()-tic

# Check Bayes estimate
colMeans(out)

# Plot trace plots
matplot(out,typ="l",lty=1)

test_bayesprobit.f90

! Written by Dries F. Benoit 
! Faculty of economics and business administration
! Ghent University - BELGIUM

! Draws from the posterior distribution of the Bayesian probit model
! with normal prior on the regression coefficients.

! Algorithm based on:
! Albert and Chib (1993). Bayesian Analysis of Binary and Polychotomous
! Response Data. Journal of the American Statistical Association, 
! 88(422), 669-679.

! Licence GLPv3 (or higher)

program test_bayesprobit

        implicit none

        ! Precision statement
        integer, parameter :: dp = kind(1.0d0)

        integer :: i
        integer, dimension(5000) :: y
        real(dp) :: tic, toc 
        real(dp), dimension(3) :: b0, bayesest 
        real(dp), dimension(3,3) :: bb0 
        real(dp), dimension(5000,3) :: x, betadraw

        ! Read in the data
        open(10, file="probitdata.csv",access='sequential',form="formatted")
        read(10,*) ! Skip first row

        do i = 1,5000
                read(10,*) y(i), x(i,:)
        enddo

        ! Prior
        b0 = 0.0_dp
        bb0 = 0.0_dp
        do i = 1,3
                bb0(i,i) = 1000.0_dp
        enddo

        ! Execute algorithm
        call cpu_time(tic)
        call binprobbayes(y,x,b0,bb0,5000,1,5000,3,betadraw)
        call cpu_time(toc)

        ! Write time
        print *, 'Time in seconds:',toc-tic

        ! Write result 
        do i = 1,3
                bayesest(i) = sum(betadraw(:,i))/5000.0_dp
        enddo
        print *, 'Bayes estimate:', bayesest
end

test_bayesprobit.jl

# Source the algorithms
include("binprobbayes.jl")
include("rtnorm.jl")

# Read in the data
using DataFrames
df = readtable("probitdata.csv")

# Execute algorithm
tic();
out = binprobbayes(convert(Array,df[:y]),convert(Array,df[collect(2:4)]),[0.; 0.; 0.],eye(3)*1000.,5000,1);
toc();

# Check Bayes estimate
mean(out,1)

# Plot trace plots
using PyPlot
plot(out)