# Code for simulations to produce Table 1 in Humphreys and Laver 2009. # Written for R 2.9 # Table 1 library(MASS) # MAKE A VARIANCE COVARIANCE MATRIX ################################################################# vcv1 = function(dims, rho){ V=matrix(rho, dims, dims) V[row(V)==col(V)] <- 1 V} vcv1(3,0) # THIS FUNCTION DRAWS A SET OF NORMALIZED IDEALS # ROWS ARE PEOPLE< COLS DIMENSIONS ################################################################# draw <- function(dims, n, vcv){ b1<- t(mvrnorm(n, numeric(dims), vcv)) t(sapply(1:dims, function(i,b) b[i,] - median(b[i,]), b=b1)) } #example draw(4,3,vcv1(4,0)) # valuation of a point ################################################################# v <- function(x,point,n){ sapply(1:n, function(i) colSums(abs(x[,i] - point))) } #example x= draw(4,3,vcv1(4,0)) x v(x, matrix(0,4), 3) # majority preference for a point ################################################################# m <- function(x,point,n){ (sum((v(x,point,n) < colSums(abs(x)))) )>n/2 } # simulation for one preference profile against "small" moves ################################################################# beatable <- function(x, sims, n, dims){ i=1 y=0 while (i<=sims){ y=y+m(x, matrix(runif(dims, -.001, .001)), n) i=i+1 + sims*y } y} # simulation across preference profiles, returns probability of being unbeatable ################################################################# share <- function(sims_pref, sims_points, n, dims, rho){ vcv = vcv1(dims, rho) z<-matrix(-100, sims_pref) for (i in 1:sims_pref){ x <- draw(dims, n, vcv) z[i]=beatable(x, sims_points, n, dims) } 1-mean(z) } # Illustration ################################################################# sims=100 n=3 dims=2 point = t(t(runif(dims, -.01, .01))) x=draw(dims,n,vcv1(dims, 0)) x point (v(x,point,n) - colSums(abs(x)))<0 m(x, point, n) beatable(x, sims, n, dims) share(100, 100, n, dims, 0) n=5 dims=3 point = t(t(runif(dims, -.01, .01))) x=draw(dims,n,vcv1(dims, 0)) x point (v(x,point,n) - colSums(abs(x)))<0 m(x, point, n) beatable(x, sims, n, dims) share(100, 100, n, dims, 0) ################################################################# # RESULTS ################################################################# s_pref = 5000 # The number of preference configurations to examine s_points = 5000 # The number of deviations to examine for each N = seq(3,15,2) D = seq(2,6,1) OUTPUT=OUTPUT2=matrix(-1, 7, 5) r = 0 for (i in 1:7){ for (j in 1:5){ OUTPUT[i, j] = share(s_pref, s_points, N[i], D[j], r) } } OUTPUT save(OUTPUT, file= "C:/Documents and Settings/Macartan/My Documents/My Dropbox/M/v/pp/CB/20091024_OUTPUT") r = .9 for (i in 1:7){ for (j in 1:5){ OUTPUT2[i, j] = share(s_pref, s_points, N[i], D[j], r) } } OUTPUT2 save(OUTPUT2, file= "C:/Documents and Settings/Macartan/My Documents/My Dropbox/M/v/pp/CB/20091024_OUTPUT2")