First Commit - changed from 'cgen', fresh git history

This commit is contained in:
Brian Albert Monroe 2017-07-11 17:02:30 +02:00
commit e5a800c822
18 changed files with 1649 additions and 0 deletions

2
.Rbuildignore Normal file
View File

@ -0,0 +1,2 @@
^.*\.Rproj$
^\.Rproj\.user$

4
.gitignore vendored Normal file
View File

@ -0,0 +1,4 @@
.Rproj.user
.Rhistory
.RData
tags

16
DESCRIPTION Normal file
View File

@ -0,0 +1,16 @@
Package: rcgen
Title: Generate Chocies for Popular Experimental Economics Instruments
Version: 0.2.0.0000
Authors@R: person("Brian", "Monroe", email = "brianalbertmonroe@gmail.com", role = c("aut", "cre"))
Description: Generates chocies from popular experimental economics instruments. This package makes simulated data that can then be used to test various estimation techniques.
Depends:
R (>= 3.3.1),
RcppArmadillo
License: MIT
Encoding: UTF-8
LazyData: true
LinkingTo: Rcpp, RcppArmadillo
Imports:
Rcpp,
RcppArmadillo
RoxygenNote: 6.0.1

6
NAMESPACE Normal file
View File

@ -0,0 +1,6 @@
# Generated by roxygen2: do not edit by hand
export(genChoice)
import(RcppArmadillo)
importFrom(Rcpp,sourceCpp)
useDynLib(rcgen)

11
R/RcppExports.R Normal file
View File

@ -0,0 +1,11 @@
# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
genEUTcpp <- function(par, dempars, N, A, B, pA, pB, Min, Max, ufunc, lo, unif, choice_dbug) {
.Call('rcgen_genEUTcpp', PACKAGE = 'rcgen', par, dempars, N, A, B, pA, pB, Min, Max, ufunc, lo, unif, choice_dbug)
}
genRDU <- function(par, dempars, N, A, B, pA, pB, Min, Max, type, ufunc, unif, lo, choice_dbug) {
.Call('rcgen_genRDU', PACKAGE = 'rcgen', par, dempars, N, A, B, pA, pB, Min, Max, type, ufunc, unif, lo, choice_dbug)
}

4
R/extra.R Normal file
View File

@ -0,0 +1,4 @@
#' @useDynLib rcgen
#' @importFrom Rcpp sourceCpp
#' @import RcppArmadillo
NULL

171
R/genfuns.r Normal file
View File

@ -0,0 +1,171 @@
CRRA <- function(x,r){
ifelse( r == 1, log(x), (x^(1-r)) / (1-r) )
}
context <- function(x, Max){
# All Options need to have equal number of outcomes, and we're only handling binary lottery pairs here
# So the first half of the x vector will be probabilities and the second half outcomes.
onum <- length(x) / 2
base <- ifelse(Max, -9999999999, 999999999)
y <- ifelse(x[1:onum]>0, x[(onum+1):length(x)] , base)
ifelse(Max, max(y), min(y))
}
# Supply the standard deviations and the correlation coefficients and get back
# the covariance matrix
mkcov <- function(sd,rho){
z <- .5 * ((8*length(rho) + 1)^.5 +1)
RHO <- diag(x=1, nrow=z, ncol=z)
count <- 0
for(row in 1:(z-1)){
for(col in (row+1):z){
count <- count + 1
RHO[row,col] <- rho[count]
}
}
RHO[lower.tri(RHO)] <- RHO[upper.tri(RHO)]
sigma <- diag(x=1, nrow=z, ncol=z)
# Fill out the covariance matrix
for(row in 1:z){
for(col in 1:z){
sigma[row,col] = sd[row] * sd[col] * RHO[row,col];
}
}
sigma
}
genEUT <- function(par, dempars, N, A, B, pA, pB, Min, Max){
rm <- par[1] # CRRA Mean
rs <- par[2] # CRRA Standard Deviation
um <- par[3] # Fechner Mean
us <- par[4] # Fechner Standard Deviation
rh <- par[5] # Rho Mean
# Fechner will use a gamma distribution, so need to back out shape and
# scale parameters
means <- c(rm,um)
sds <- c(rs,us)
sigma <- mkcov(sds,rh)
dists <- MASS::mvrnorm(n=N, mu=means, Sigma=sigma)
p1 <- pnorm(dists[,1], mean=means[1], sd=sds[1])
p2 <- pnorm(dists[,2], mean=means[2], sd=sds[2])
pvals <- cbind(p1,p2)
k <- (um^2)/(us^2)
t <- (us^2)/um
r <- qnorm(pvals[,1], mean=rm, sd=rs)
mu <- qgamma(pvals[,2], shape=k, scale=t)
d <- data.frame(choice=rep(0, nrow(A)))
for(i in 0:(ncol(A)-1)){
d[[paste0("A",i)]] <- A[,i+1]
}
for(i in 0:(ncol(A)-1)){
d[[paste0("pA",i)]] <- pA[,i+1]
}
for(i in 0:(ncol(B)-1)){
d[[paste0("B",i)]] <- B[,i+1]
}
for(i in 0:(ncol(B)-1)){
d[[paste0("pB",i)]] <- pB[,i+1]
}
d$Min <- Min
d$Max <- Max
e <- d
e$ID <- 1
rval <- r[1]
uval <- mu[1]
while( is.na(rval) | is.infinite(rval)){
rval <- rnorm(1,mean = rm, sd = rs)
}
while( is.na(uval) | is.infinite(uval)){
uval <- rgamma(1,shape = k, scale = t)
}
e$r <- rval
e$mu <- uval
D <- e
for( i in 2:N){
e <- d
e$ID <- i
# grab one value from the distribution
rval <- r[i]
uval <- mu[i]
while( is.na(rval) | is.infinite(rval)){
rval <- rnorm(1,mean = rm, sd = rs)
}
while( is.na(uval) | is.infinite(uval)){
uval <- rgamma(1,shape = k, scale = t)
}
e$r <- rval
e$mu <- uval
D <- rbind(D,e)
}
ctx <- CRRA(D$Max,D$r) - CRRA(D$Min,D$r)
UA <- 0
UB <- 0
for(i in 1:ncol(A)){
UA <- UA + (pA[,i] * CRRA(A[,i], D$r))
UB <- UB + (pB[,i] * CRRA(B[,i], D$r))
}
UB.1 <- (UB/ctx/D$mu) - (UA/ctx/D$mu)
UA.1 <- 0
# Have things gone haywire because of the computer's inability to handle
# numbers bigger than ~3e310 ?
c.N <- is.nan(UB.1) | is.infinite(UB.1)
Aprob <- ifelse( c.N , # Are we dealing with an insane number?
# yes:
ifelse( UB > UA , 0 , 1 ) ,
# no, but are we making an insane number via exp?
ifelse( UB.1 > 709 , 0 ,
ifelse( UB.1 < -709, 1 , exp(UA.1) / (exp(UA.1) + exp(UB.1)) )
)
)
Bprob <- 1 - Aprob
# Random uniform number
rand <- runif(nrow(D))
# This is a great R function, ifelse collapses would would otherwise
# potentially be several lines of code into a very readable one line
# statement.
D$choice <- ifelse(Aprob > rand, 0, 1)
return(D)
}
genEUT <- compiler::cmpfun(genEUT)

336
R/geninst.r Normal file
View File

@ -0,0 +1,336 @@
genHL <- function(par, N, dempars=list()) {
# Holt & Laury (2002)
# 10 lotteries
A0 <- rep(1.60, 10)
A1 <- rep(2.00, 10)
B0 <- rep(0.10, 10)
B1 <- rep(3.85, 10)
pA0 <- seq(from=0.9, to=0, by= -.1)
pB0 <- seq(from=0.9, to=0, by= -.1)
pA1 <- seq(from=0.1, to=1, by= .1)
pB1 <- seq(from=0.1, to=1, by= .1)
A <- cbind(A0,A1)
B <- cbind(B0,B1)
pA <- cbind(pA0, pA1)
pB <- cbind(pB0, pB1)
Max <- apply( cbind(pA,pB,A,B) ,1, context, Max = T )
Min <- apply( cbind(pA,pB,A,B) ,1, context, Max = F )
type <- rep("HL", length(A0))
list(par = par, dempars = dempars, N = N, A = A, B = B, pA = pA, pB = pB, Min = Min, Max = Max, inst_type = type)
}
genHO <- function(par, N, dempars=list()) {
# Hey & Orme (1994) - Table I
# 100 Lotteries
A0 <- c(rep(0 , 25), rep(0 , 25), rep(0 , 25), rep(10, 25))
A1 <- c(rep(10, 25), rep(10, 25), rep(20, 25), rep(20, 25))
A2 <- c(rep(20, 25), rep(30, 25), rep(30, 25), rep(30, 25))
B0 <- c(rep(0 , 25), rep(0 , 25), rep(0 , 25), rep(10, 25))
B1 <- c(rep(10, 25), rep(10, 25), rep(20, 25), rep(20, 25))
B2 <- c(rep(20, 25), rep(30, 25), rep(30, 25), rep(30, 25))
pA0 <- c(.625, .375, .000, .125, .500, .250, .250, .250, .125, .125, .125, .250, .625, .125, .125, .375, .000, .500, .750, .250, .000, .000, .250, .500, .250)
pA1 <- c(.000, .625, 1.00, .750, .375, .750, .625, .250, .375, .250, .875, .750, .375, .500, .750, .375, .750, .125, .000, .375, .875, .625, .500, .500, .500)
pA2 <- c(.375, .000, .000, .125, .125, .000, .125, .500, .500, .625, .000, .000, .000, .375, .125, .250, .250, .375, .250, .375, .125, .375, .250, .000, .250)
pA0 <- rep(pA0, 4)
pA1 <- rep(pA1, 4)
pA2 <- rep(pA2, 4)
pB0 <- c(.375, .500, .125, .250, .625, .375, .375, .125, .000, .000, .250, .500, .750, .250, .375, .500, .125, .375, .625, .375, .125, .125, .125, .625, .375)
pB1 <- c(.625, .250, .500, .500, .125, .000, .250, .625, 1.00, .500, .625, .000, .125, .000, .125, .125, .375, .500, .375, .000, .625, .250, .875, .125, .250)
pB2 <- c(.000, .250, .375, .250, .250, .625, .375, .250, .000, .500, .125, .500, .125, .750, .500, .375, .500, .125, .000, .625, .250, .625, .000, .250, .375)
pB0 <- rep(pB0, 4)
pB1 <- rep(pB1, 4)
pB2 <- rep(pB2, 4)
A <- cbind(A0, A1, A2)
B <- cbind(B0, B1, B2)
pA <- cbind(pA0, pA1, pA2)
pB <- cbind(pB0, pB1, pB2)
Max <- apply(cbind(pA, pB, A, B), 1, context, Max = T )
Min <- apply(cbind(pA, pB, A, B), 1, context, Max = F )
type <- rep("HO", length(A0))
list(par = par, dempars = dempars, N = N, A = A, B = B, pA = pA, pB = pB, Min = Min, Max = Max, inst_type = type)
}
genHNG <- function(par, N, dempars=list()) {
# Harrison & NG (2016)
# 80 Lotteries
pA0 <- c(0,0.15,0,0.15,0.15,0.6,0.6,0.9,0,0.1,0.5,0,0.5,0.4,0.4,0.9,0,0.1,0,0,0.7,0.7,0.6,0.75,0.1,0.2,0.4,0.1,0.5,0.6,0.4,0.8,0,0.1,0,0,0.5,0.55,0.5,0.7,0.03,0.18,0.27,0.12,0.06,0.54,0.18,0.84,0.08,0.2,0.1,0.08,0.65,0.65,0.48,0.83,0.08,0.38,0.38,0.02,0.62,0.66,0.52,0.81,0.08,0.25,0.35,0.18,0.53,0.55,0.48,0.78,0.08,0.22,0.2,0.02,0.48,0.45,0.54,0.65)
pA1 <- c(0.25,0.25,0.5,0.25,0.75,0,0,0,0.2,0.8,0,1,0.4,0.6,0.6,0,0.25,0.75,1,1,0,0,0.25,0.25,0,0.6,0,0.9,0.3,0,0.6,0,0.4,0.6,1,1,0.2,0,0.2,0,0.2,0.2,0.05,0.3,0.9,0.1,0.7,0.1,0.04,0.6,0.8,0.84,0.1,0.1,0.44,0.14,0.05,0.05,0.05,0.95,0.2,0.1,0.45,0.1,0.06,0.45,0.15,0.66,0.21,0.15,0.36,0.06,0.08,0.12,0.2,0.92,0.28,0.4,0.04,0.2)
pA2 <- c(0.75,0.6,0.5,0.6,0.1,0.4,0.4,0.1,0.8,0.1,0.5,0,0.1,0,0,0.1,0.75,0.15,0,0,0.3,0.3,0.15,0,0.9,0.2,0.6,0,0.2,0.4,0,0.2,0.6,0.3,0,0,0.3,0.45,0.3,0.3,0.77,0.62,0.68,0.58,0.04,0.36,0.12,0.06,0.88,0.2,0.1,0.08,0.25,0.25,0.08,0.03,0.87,0.57,0.57,0.03,0.18,0.24,0.03,0.09,0.86,0.3,0.5,0.16,0.26,0.3,0.16,0.16,0.84,0.66,0.6,0.06,0.24,0.15,0.42,0.15)
pB0 <- c(0.15,0.3,0.3,0,0,0,0.15,0.75,0.1,0.5,0,0.1,0.7,0.7,0.5,0.8,0.1,0.4,0.4,0.1,0.6,0.5,0.5,0.85,0,0.4,0.1,0.2,0.6,0.4,0.5,0.7,0.1,0.25,0.25,0.1,0.4,0.4,0.55,0.6,0.12,0.27,0.03,0.03,0.12,0.06,0.54,0.78,0.05,0.45,0.45,0.05,0.55,0.45,0.44,0.88,0.04,0.14,0.04,0.08,0.68,0.54,0.58,0.77,0.02,0.35,0.15,0.12,0.58,0.45,0.42,0.72,0.02,0.13,0.1,0.08,0.44,0.5,0.52,0.68)
pB1 <- c(0,0,0,0.5,1,1,0.75,0.25,0,0,1,0.8,0,0,0.4,0.2,0,0,0,0.75,0.25,0.5,0.5,0,0.3,0,0.9,0.6,0,0.6,0.3,0.3,0,0,0,0.6,0.6,0.6,0,0.4,0.05,0.05,0.45,0.45,0.8,0.9,0.1,0.2,0.1,0.1,0.1,0.9,0.3,0.5,0.52,0.04,0.15,0.65,0.9,0.8,0.05,0.4,0.3,0.2,0.24,0.15,0.75,0.84,0.06,0.45,0.54,0.24,0.32,0.48,0.6,0.68,0.44,0.2,0.12,0.08)
pB2 <- c(0.85,0.7,0.7,0.5,0,0,0.1,0,0.9,0.5,0,0.1,0.3,0.3,0.1,0,0.9,0.6,0.6,0.15,0.15,0,0,0.15,0.7,0.6,0,0.2,0.4,0,0.2,0,0.9,0.75,0.75,0.3,0,0,0.45,0,0.83,0.68,0.52,0.52,0.08,0.04,0.36,0.02,0.85,0.45,0.45,0.05,0.15,0.05,0.04,0.08,0.81,0.21,0.06,0.12,0.27,0.06,0.12,0.03,0.74,0.5,0.1,0.04,0.36,0.1,0.04,0.04,0.66,0.39,0.3,0.24,0.12,0.3,0.36,0.24)
A0 <- rep(10, 80)
A1 <- rep(30, 80)
A2 <- rep(50, 80)
B0 <- rep(10, 80)
B1 <- rep(30, 80)
B2 <- rep(50, 80)
A <- cbind(A0,A1,A2)
B <- cbind(B0,B1,B2)
pA <- cbind(pA0, pA1, pA2)
pB <- cbind(pB0, pB1, pB2)
Max <- apply( cbind(pA,pB,A,B) ,1, context, Max = T )
Min <- apply( cbind(pA,pB,A,B) ,1, context, Max = F )
type <- rep("HNG", length(A0))
list(par = par, dempars = dempars, N = N, A = A, B = B, pA = pA, pB = pB, Min = Min, Max = Max, inst_type = type)
}
genHNG.ins <- function(par, N, dempars=list()) {
# Harrison & NG (2016)
# 24 lotteries
lotnum <- 24
B0 <- c(19.8, 19.6, 19.4, 19.2, 19, 18.8, 18.6, 18.4, 18.2, 18, 17.8, 17.6, 17.4, 17.2, 17, 16.8, 16.6, 16.4, 16.2, 16, 15.8, 15.6, 15.4, 15.2)
B1 <- rep(1, lotnum)
B2 <- rep(1, lotnum)
A0 <- rep(5, lotnum)
A1 <- rep(20, lotnum)
A2 <- rep(1, lotnum)
pA0 <- rep(0.1, lotnum)
pA1 <- rep(0.9, lotnum)
pA2 <- rep(0, lotnum)
pB0 <- rep(1, lotnum)
pB1 <- rep(0, lotnum)
pB2 <- rep(0, lotnum)
A <- cbind(A0, A1, A2)
B <- cbind(B0, B1, B2)
pA <- cbind(pA0, pA1, pA2)
pB <- cbind(pB0, pB1, pB2)
Max <- apply( cbind(pA,pB,A,B) ,1, context, Max = T )
Min <- apply( cbind(pA,pB,A,B) ,1, context, Max = F )
type <- rep("HNG.ins", length(A0))
list(par = par, dempars = dempars, N = N, A = A, B = B, pA = pA, pB = pB, Min = Min, Max = Max, inst_type = type)
}
genHNG.lo <- function(par, N, dempars=list()) {
# Harrison & NG (2016)
# 80 lotteries
pA0 <- c(0,0.15,0,0.15,0.15,0.6,0.6,0.9,0,0.1,0.5,0,0.5,0.4,0.4,0.9,0,0.1,0,0,0.7,0.7,0.6,0.75,0.1,0.2,0.4,0.1,0.5,0.6,0.4,0.8,0,0.1,0,0,0.5,0.55,0.5,0.7,0.03,0.18,0.27,0.12,0.06,0.54,0.18,0.84,0.08,0.2,0.1,0.08,0.65,0.65,0.48,0.83,0.08,0.38,0.38,0.02,0.62,0.66,0.52,0.81,0.08,0.25,0.35,0.18,0.53,0.55,0.48,0.78,0.08,0.22,0.2,0.02,0.48,0.45,0.54,0.65)
pA1 <- c(0.25,0.25,0.5,0.25,0.75,0,0,0,0.2,0.8,0,1,0.4,0.6,0.6,0,0.25,0.75,1,1,0,0,0.25,0.25,0,0.6,0,0.9,0.3,0,0.6,0,0.4,0.6,1,1,0.2,0,0.2,0,0.2,0.2,0.05,0.3,0.9,0.1,0.7,0.1,0.04,0.6,0.8,0.84,0.1,0.1,0.44,0.14,0.05,0.05,0.05,0.95,0.2,0.1,0.45,0.1,0.06,0.45,0.15,0.66,0.21,0.15,0.36,0.06,0.08,0.12,0.2,0.92,0.28,0.4,0.04,0.2)
pA2 <- c(0.75,0.6,0.5,0.6,0.1,0.4,0.4,0.1,0.8,0.1,0.5,0,0.1,0,0,0.1,0.75,0.15,0,0,0.3,0.3,0.15,0,0.9,0.2,0.6,0,0.2,0.4,0,0.2,0.6,0.3,0,0,0.3,0.45,0.3,0.3,0.77,0.62,0.68,0.58,0.04,0.36,0.12,0.06,0.88,0.2,0.1,0.08,0.25,0.25,0.08,0.03,0.87,0.57,0.57,0.03,0.18,0.24,0.03,0.09,0.86,0.3,0.5,0.16,0.26,0.3,0.16,0.16,0.84,0.66,0.6,0.06,0.24,0.15,0.42,0.15)
pB0 <- c(0.15,0.3,0.3,0,0,0,0.15,0.75,0.1,0.5,0,0.1,0.7,0.7,0.5,0.8,0.1,0.4,0.4,0.1,0.6,0.5,0.5,0.85,0,0.4,0.1,0.2,0.6,0.4,0.5,0.7,0.1,0.25,0.25,0.1,0.4,0.4,0.55,0.6,0.12,0.27,0.03,0.03,0.12,0.06,0.54,0.78,0.05,0.45,0.45,0.05,0.55,0.45,0.44,0.88,0.04,0.14,0.04,0.08,0.68,0.54,0.58,0.77,0.02,0.35,0.15,0.12,0.58,0.45,0.42,0.72,0.02,0.13,0.1,0.08,0.44,0.5,0.52,0.68)
pB1 <- c(0,0,0,0.5,1,1,0.75,0.25,0,0,1,0.8,0,0,0.4,0.2,0,0,0,0.75,0.25,0.5,0.5,0,0.3,0,0.9,0.6,0,0.6,0.3,0.3,0,0,0,0.6,0.6,0.6,0,0.4,0.05,0.05,0.45,0.45,0.8,0.9,0.1,0.2,0.1,0.1,0.1,0.9,0.3,0.5,0.52,0.04,0.15,0.65,0.9,0.8,0.05,0.4,0.3,0.2,0.24,0.15,0.75,0.84,0.06,0.45,0.54,0.24,0.32,0.48,0.6,0.68,0.44,0.2,0.12,0.08)
pB2 <- c(0.85,0.7,0.7,0.5,0,0,0.1,0,0.9,0.5,0,0.1,0.3,0.3,0.1,0,0.9,0.6,0.6,0.15,0.15,0,0,0.15,0.7,0.6,0,0.2,0.4,0,0.2,0,0.9,0.75,0.75,0.3,0,0,0.45,0,0.83,0.68,0.52,0.52,0.08,0.04,0.36,0.02,0.85,0.45,0.45,0.05,0.15,0.05,0.04,0.08,0.81,0.21,0.06,0.12,0.27,0.06,0.12,0.03,0.74,0.5,0.1,0.04,0.36,0.1,0.04,0.04,0.66,0.39,0.3,0.24,0.12,0.3,0.36,0.24)
A0 <- rep(00, 80)
A1 <- rep(20, 80)
A2 <- rep(40, 80)
B0 <- rep(00, 80)
B1 <- rep(20, 80)
B2 <- rep(40, 80)
A <- cbind(A0,A1,A2)
B <- cbind(B0,B1,B2)
pA <- cbind(pA0, pA1, pA2)
pB <- cbind(pB0, pB1, pB2)
Max <- apply( cbind(pA,pB,A,B) ,1, context, Max = T )
Min <- apply( cbind(pA,pB,A,B) ,1, context, Max = F )
type <- rep("HNG.lot", length(A0))
list(par = par, dempars = dempars, N = N, A = A, B = B, pA = pA, pB = pB, Min = Min, Max = Max, inst_type = type)
}
genLMS20 <- function(par, N, dempars=list()) {
# Loomes, Moffat & Sugden (2002)
# 90 lotteries
pA0 <- c(0.15, 0.30, 0.30, 0.15, 0.15, 0.60, 0.60, 0.90, 0.10, 0.50, 0.50, 0.10, 0.70, 0.70, 0.50, 0.90, 0.10, 0.40, 0.40, 0.10, 0.70, 0.70, 0.60, 0.85, 0.10, 0.40, 0.40, 0.20, 0.60, 0.60, 0.50, 0.80, 0.10, 0.20, 0.20, 0.10, 0.35, 0.40, 0.40, 0.60, 0.00, 0.55, 0.80, 0.10, 0.70)
pA0 <- c(pA0, pA0)
pA1 <- c(0.00, 0.00, 0.00, 0.25, 0.75, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.80, 0.00, 0.00, 0.40, 0.00, 0.00, 0.00, 0.00, 0.75, 0.00, 0.00, 0.25, 0.00, 0.00, 0.00, 0.00, 0.60, 0.00, 0.00, 0.30, 0.00, 0.00, 0.00, 0.00, 0.50, 0.25, 0.00, 0.00, 0.00, 0.25, 0.20, 0.00, 0.75, 0.30)
pA1 <- c(pA1, pA1)
pA2 <- c(0.85, 0.70, 0.70, 0.60, 0.10, 0.40, 0.40, 0.10, 0.90, 0.50, 0.50, 0.10, 0.30, 0.30, 0.10, 0.10, 0.90, 0.60, 0.60, 0.15, 0.30, 0.30, 0.15, 0.15, 0.90, 0.60, 0.60, 0.20, 0.40, 0.40, 0.20, 0.20, 0.90, 0.80, 0.80, 0.40, 0.40, 0.60, 0.60, 0.40, 0.75, 0.25, 0.20, 0.15, 0.00)
pA2 <- c(pA2, pA2)
pB0 <- c(0.00, 0.15, 0.00, 0.00, 0.00, 0.00, 0.15, 0.75, 0.00, 0.10, 0.00, 0.00, 0.50, 0.40, 0.40, 0.80, 0.00, 0.10, 0.00, 0.00, 0.60, 0.50, 0.50, 0.75, 0.00, 0.20, 0.10, 0.10, 0.50, 0.40, 0.40, 0.70, 0.00, 0.10, 0.00, 0.00, 0.25, 0.25, 0.35, 0.50, 0.00, 0.65, 0.85, 0.15, 0.75)
pB0 <- c(pB0, pB0)
pB1 <- c(0.25, 0.25, 0.50, 0.50, 1.00, 1.00, 0.75, 0.25, 0.20, 0.80, 1.00, 1.00, 0.40, 0.60, 0.60, 0.20, 0.25, 0.75, 1.00, 1.00, 0.25, 0.50, 0.50, 0.25, 0.30, 0.60, 0.90, 0.90, 0.30, 0.60, 0.60, 0.30, 0.50, 0.50, 1.00, 1.00, 0.75, 0.75, 0.25, 0.50, 0.30, 0.15, 0.00, 0.75, 0.25)
pB1 <- c(pB1, pB1)
pB2 <- c(0.75, 0.60, 0.50, 0.50, 0.00, 0.00, 0.10, 0.00, 0.80, 0.10, 0.00, 0.00, 0.10, 0.00, 0.00, 0.00, 0.75, 0.15, 0.00, 0.00, 0.15, 0.00, 0.00, 0.00, 0.70, 0.20, 0.00, 0.00, 0.20, 0.00, 0.00, 0.00, 0.50, 0.40, 0.00, 0.00, 0.00, 0.00, 0.40, 0.00, 0.70, 0.20, 0.15, 0.10, 0.00)
pB2 <- c(pB2, pB2)
A0 <- rep(00, 90)
A1 <- rep(10, 90)
A2 <- rep(20, 90)
B0 <- rep(00, 90)
B1 <- rep(10, 90)
B2 <- rep(20, 90)
A <- cbind(A0, A1, A2)
B <- cbind(B0, B1, B2)
pA <- cbind(pA0, pA1, pA2)
pB <- cbind(pB0, pB1, pB2)
Max <- apply( cbind(pA,pB,A,B) ,1, context, Max = T )
Min <- apply( cbind(pA,pB,A,B) ,1, context, Max = F )
type <- rep("LMS20", length(A0))
list(par = par, dempars = dempars, N = N, A = A, B = B, pA = pA, pB = pB, Min = Min, Max = Max, inst_type = type)
}
genLMS30 <- function(par, N, dempars=list()) {
# Loomes, Moffat & Sugden (2002)
# 90 lotteries
pA0 <- c(0.15, 0.30, 0.30, 0.15, 0.15, 0.60, 0.60, 0.90, 0.10, 0.50, 0.50, 0.10, 0.70, 0.70, 0.50, 0.90, 0.10, 0.40, 0.40, 0.10, 0.70, 0.70, 0.60, 0.85, 0.10, 0.40, 0.40, 0.20, 0.60, 0.60, 0.50, 0.80, 0.10, 0.20, 0.20, 0.10, 0.35, 0.10, 0.25, 0.25, 0.10, 0.50, 0.55, 0.55, 0.70)
pA0 <- c(pA0, pA0)
pA1 <- c(0.00, 0.00, 0.00, 0.25, 0.75, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.80, 0.00, 0.00, 0.40, 0.00, 0.00, 0.00, 0.00, 0.75, 0.00, 0.00, 0.25, 0.00, 0.00, 0.00, 0.00, 0.60, 0.00, 0.00, 0.30, 0.00, 0.00, 0.00, 0.00, 0.50, 0.25, 0.00, 0.00, 0.00, 0.60, 0.20, 0.00, 0.00, 0.00)
pA1 <- c(pA1, pA1)
pA2 <- c(0.85, 0.70, 0.70, 0.60, 0.10, 0.40, 0.40, 0.10, 0.90, 0.50, 0.50, 0.10, 0.30, 0.30, 0.10, 0.10, 0.90, 0.60, 0.60, 0.15, 0.30, 0.30, 0.15, 0.15, 0.90, 0.60, 0.60, 0.20, 0.40, 0.40, 0.20, 0.20, 0.90, 0.80, 0.80, 0.40, 0.40, 0.90, 0.75, 0.75, 0.30, 0.30, 0.45, 0.45, 0.30)
pA2 <- c(pA2, pA2)
pB0 <- c(0.00, 0.15, 0.00, 0.00, 0.00, 0.00, 0.15, 0.75, 0.00, 0.10, 0.00, 0.00, 0.50, 0.40, 0.40, 0.80, 0.00, 0.10, 0.00, 0.00, 0.60, 0.50, 0.50, 0.75, 0.00, 0.20, 0.10, 0.10, 0.50, 0.40, 0.40, 0.70, 0.00, 0.10, 0.00, 0.00, 0.25, 0.00, 0.10, 0.00, 0.00, 0.40, 0.40, 0.50, 0.60)
pB0 <- c(pB0, pB0)
pB1 <- c(0.25, 0.25, 0.50, 0.50, 1.00, 1.00, 0.75, 0.25, 0.20, 0.80, 1.00, 1.00, 0.40, 0.60, 0.60, 0.20, 0.25, 0.75, 1.00, 1.00, 0.25, 0.50, 0.50, 0.25, 0.30, 0.60, 0.90, 0.90, 0.30, 0.60, 0.60, 0.30, 0.50, 0.50, 1.00, 1.00, 0.75, 0.40, 0.60, 1.00, 1.00, 0.60, 0.60, 0.20, 0.40)
pB1 <- c(pB1, pB1)
pB2 <- c(0.75, 0.60, 0.50, 0.50, 0.00, 0.00, 0.10, 0.00, 0.80, 0.10, 0.00, 0.00, 0.10, 0.00, 0.00, 0.00, 0.75, 0.15, 0.00, 0.00, 0.15, 0.00, 0.00, 0.00, 0.70, 0.20, 0.00, 0.00, 0.20, 0.00, 0.00, 0.00, 0.50, 0.40, 0.00, 0.00, 0.00, 0.60, 0.30, 0.00, 0.00, 0.00, 0.00, 0.30, 0.00)
pB2 <- c(pB2, pB2)
A0 <- rep(00, 90)
A1 <- rep(10, 90)
A2 <- rep(30, 90)
B0 <- rep(00, 90)
B1 <- rep(10, 90)
B2 <- rep(30, 90)
A <- cbind(A0, A1, A2)
B <- cbind(B0, B1, B2)
pA <- cbind(pA0, pA1, pA2)
pB <- cbind(pB0, pB1, pB2)
Max <- apply( cbind(pA,pB,A,B) ,1, context, Max = T )
Min <- apply( cbind(pA,pB,A,B) ,1, context, Max = F )
type <- rep("LMS30", length(A0))
list(par = par, dempars = dempars, N = N, A = A, B = B, pA = pA, pB = pB, Min = Min, Max = Max, inst_type = type)
}
genSH <- function(par, N, dempars=list()) {
#Schmidt & Hey (2004)
# 150 Lotteries
pA0 <- c(.000, .750, .300, .000, .000, .000, .500, .000, .800, .200, .000, .500, .000, .000, .200, .100, .000, .500, .000, .000, .000, .500, .250, .000, .500, .000, .000, .250)
pA1 <- c(.000, .000, .600, .600, 1.00, .500, .500, .000, .000, .000, .200, .100, .200, .100, .800, .400, .400, .200, .200, .200, .000, .000, .500, .500, .250, .250, .000, .250)
pA2 <- c(1.00, .250, .100, .100, .000, .500, .000, .700, .140, .740, .800, .400, .600, .300, .000, .500, .600, .300, .300, .700, .500, .500, .250, .000, .000, .500, .750, .500)
pA3 <- c(.000, .000, .000, .300, .000, .000, .000, .300, .060, .060, .000, .000, .200, .600, .000, .000, .000, .000, .500, .100, .500, .000, .000, .500, .250, .250, .250, .000)
pA0 <- c(pA0, pA0, pA0, pA0, pA0)
pA1 <- c(pA1, pA1, pA1, pA1, pA1)
pA2 <- c(pA2, pA2, pA2, pA2, pA2)
pA3 <- c(pA3, pA3, pA3, pA3, pA3)
pB0 <- c(.200, .800, .320, .020, .700, .350, .850, .150, .830, .230, .000, .500, .200, .100, .800, .400, .400, .700, .200, .200, .100, .600, .300, .200, .600, .000, .000, .250)
pB1 <- c(.000, .000, .600, .600, .000, .000, .000, .000, .000, .000, .500, .250, .000, .000, .000, .000, .000, .000, .000, .000, .000, .000, .500, .200, .100, .350, .100, .350)
pB2 <- c(.000, .000, .000, .000, .000, .500, .000, .000, .000, .600, .000, .000, .400, .200, .000, .500, .000, .000, .000, .400, .000, .000, .000, .000, .000, .000, .250, .000)
pB3 <- c(.800, .200, .080, .380, .300, .150, .150, .850, .170, .170, .500, .250, .400, .700, .200, .100, .600, .300, .800, .400, .900, .400, .200, .600, .300, .650, .650, .400)
pB0 <- c(pB0, pB0, pB0, pB0, pB0)
pB1 <- c(pB1, pB1, pB1, pB1, pB1)
pB2 <- c(pB2, pB2, pB2, pB2, pB2)
pB3 <- c(pB3, pB3, pB3, pB3, pB3)
A0 <- rep(00, 140)
A1 <- rep(10, 140)
A2 <- rep(30, 140)
A3 <- rep(40, 140)
B0 <- rep(00, 140)
B1 <- rep(10, 140)
B2 <- rep(30, 140)
B3 <- rep(40, 140)
A <- cbind(A0, A1, A2, A3)
B <- cbind(B0, B1, B2, B3)
pA <- cbind(pA0, pA1, pA2, pA3)
pB <- cbind(pB0, pB1, pB2, pB3)
Max <- apply( cbind(pA,pB,A,B) ,1, context, Max = T )
Min <- apply( cbind(pA,pB,A,B) ,1, context, Max = F )
type <- rep("SH", nrow(A))
list(par = par, dempars = dempars, N = N, A = A, B = B, pA = pA, pB = pB, Min = Min, Max = Max, inst_type = type)
}
genInstruments <- function(insts, par, N, dempars ) {
# We want each subject to be able to repsond to multiple instruments
# All instruments need to share the same "gen<inst>" name structure and take the same arguments
setup <- lapply(insts, function(inst, par, N, dempars){
genInst <- get(paste0("gen", inst))
genInst(par, N, dempars)
}, par = par, N = N, dempars = dempars)
# This function is useful when being called from mapply through do.call
mapping <- function(...) {
dots <- list(...)
if (is.matrix(dots[[1]])) {
cols <- lapply(dots,ncol)
mcol <- do.call(max, cols)
dots <- lapply(dots, function(x) {
xcol <- ncol(x)
if (xcol < mcol) {
namat <- matrix(NA, nrow=nrow(x), ncol = (mcol - xcol))
cbind(x,namat)
} else{
x
}
})
do.call(rbind,dots)
} else {
c(...)
}
}
# Mapply applies the the function FUN to the first element of the lists passed in the dots of mapply
# Name this FUN mapping so we can use do.call
setup$FUN <- mapping
out <- do.call(mapply, setup )
out$pA <- ifelse(is.na(out$pA), 0, out$pA)
out$pB <- ifelse(is.na(out$pB), 0, out$pB)
out$A <- ifelse(is.na(out$A), 1, out$A)
out$B <- ifelse(is.na(out$B), 1, out$B)
out$N <- out$N[1]
out$par <- out$par[1:length(par)]
out
}

118
R/main.R Normal file
View File

@ -0,0 +1,118 @@
#' Generate chocies using popular instruments from experimental economics
#'
#' @param par vector containing means and standard deviations
#' @param N Number of subjects
#' @param dempars list of covariates
#' @param inst instrument used to generate chocies: "HL", "HNG", "HNG.ins",
#' "HO", "LMS20", "LMS30"
#' @param pfunc what kind of probability weighting function, EUT, pow, invs
#' @param ufunc what kind of utility function, CRRA, pow
#' @param unif should the parameters be distributed uniformly
#' @param lo should the r-value be restricted to less than 1? (can't be used with unif)
#' @param dbug show calculated utilities and probabilities, True, False
#' @export genChoice
genChoice <- function(par, N, dempars=list(), inst = c("HL"), pfunc = "EUT", ufunc = "CRRA", unif = F, lo = F, dbug = F) {
# Different population types available
pf.available <- c("EUT", "pow", "invs", "prelec")
uf.available <- c("CRRA", "pow")
inst.available <- c("HL", "HO", "HNG", "HNG.ins", "LMS20", "LMS30", "SH", "HNG.lo")
if (! pfunc %in% pf.available) {
stop(paste("pfunc passed not available, choose from:", paste(pf.available, collapse=", ")))
}
if (! ufunc %in% uf.available) {
stop(paste("ufunc passed not available, choose from:", paste(uf.available, collapse=", ")))
}
for (i in inst) {
if (! i %in% inst.available){
stop(paste0("Instrument \"", i, "\" type passed not available, choose from: ", paste(inst.available, collapse=", ")))
}
}
# Which utility function to use
if (ufunc == "CRRA") ufunc.int = 0
else if (ufunc == "pow") ufunc.int = 1
# Change the boolian into int
dbug.int <- ifelse(dbug, 1, 0)
lo.int <- ifelse(lo, 1, 0)
unif.int <- ifelse(unif, 1, 0)
# Get the istrument stuff
setup <- genInstruments(inst, par, N, dempars)
# Reorder variables from highest to lowest for RDU, doesn't matter for EUT so do it for everyone
setup <- lapply(setup, function(x) {
if (is.matrix(x) | is.data.frame(x)) {
x[, order(colnames(x), decreasing = TRUE)]
}
else x
})
if (pfunc %in% c("pow", "invs", "prelec")) {
# Get the names of the probability variables so we can dbug with RDU showing the decision weights
pnames <- c(colnames(setup$pA), colnames(setup$pB))
pnames <- paste0("d", pnames)
} else {
pnames <- c()
}
# Add additional variables to the setup
setup$ufunc <- ufunc.int
setup$choice_dbug <- dbug.int
setup$lo <- lo.int
setup$unif <- unif.int
# Remove variables that we don't want passed to the c++ code
inst_type <- setup$inst_type
setup$inst_type <- NULL
if (pfunc == "EUT") {
if (length(par) != 5) stop("par must have 5 elements")
Dat <- do.call(genEUTcpp, setup)
pfunc.names <- c("r", "mu")
} else if (pfunc == "pow") {
if (length(par) != 9) stop("par must have 9 elements")
# There's a unified RDU function for 1 parameter RDU functions
setup$type <- 0
Dat <- do.call(genRDU, setup)
pfunc.names <- c("r", "alpha", "mu")
} else if (pfunc == "invs") {
# There's a unified RDU function for 1 parameter RDU functions
setup$type <- 1
if (length(par) != 9) stop("par must have 9 elements")
Dat <- do.call(genRDU, setup)
pfunc.names <- c("r", "alpha", "mu")
} else if (pfunc == "prelec") {
# There's a unified RDU function for 1 parameter RDU functions
setup$type <- 2
if (length(par) != 14) stop("par must have 14 elements")
Dat <- do.call(genRDU, setup)
pfunc.names <- c("r", "alpha", "beta", "mu")
}
# Setup the names of the output columns
inst.names <- c("ID", colnames(setup$A), colnames(setup$pA), colnames(setup$B), colnames(setup$pB), "Max", "Min")
if (dbug) {
inst.names <- c(inst.names, pnames, "UA", "UB", "UMin", "Umax", "CTX", "PA", "PB", "rand", "choice", pfunc.names)
} else {
inst.names <- c(inst.names, "choice", pfunc.names)
}
inst.names <- c(inst.names, "Inst", "pfunc")
# Turn Dat into a data.frame
Dat <- data.frame(Dat)
# Add the inst_type back into the Dat
Dat$Inst <- rep(inst_type, N)
# Add the pfunc
Dat$pfunc <- pfunc
# Add the names
colnames(Dat) <- inst.names
# Make the ID an integer
Dat$ID <- as.integer(Dat$ID)
Dat
}

32
man/genChoice.Rd Normal file
View File

@ -0,0 +1,32 @@
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/main.R
\name{genChoice}
\alias{genChoice}
\title{Generate chocies using popular instruments from experimental economics}
\usage{
genChoice(par, N, dempars = list(), inst = c("HL"), pfunc = "EUT",
ufunc = "CRRA", unif = F, lo = F, dbug = F)
}
\arguments{
\item{par}{vector containing means and standard deviations}
\item{N}{Number of subjects}
\item{dempars}{list of covariates}
\item{inst}{instrument used to generate chocies: "HL", "HNG", "HNG.ins",
"HO", "LMS20", "LMS30"}
\item{pfunc}{what kind of probability weighting function, EUT, pow, invs}
\item{ufunc}{what kind of utility function, CRRA, pow}
\item{unif}{should the parameters be distributed uniformly}
\item{lo}{should the r-value be restricted to less than 1? (can't be used with unif)}
\item{dbug}{show calculated utilities and probabilities, True, False}
}
\description{
Generate chocies using popular instruments from experimental economics
}

3
src/.gitignore vendored Normal file
View File

@ -0,0 +1,3 @@
*.o
*.so
*.dll

1
src/Makevars Normal file
View File

@ -0,0 +1 @@
PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)

2
src/Makevars.win Normal file
View File

@ -0,0 +1,2 @@
PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)

66
src/RcppExports.cpp Normal file
View File

@ -0,0 +1,66 @@
// Generated by using Rcpp::compileAttributes() -> do not edit by hand
// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
#include <RcppArmadillo.h>
#include <Rcpp.h>
using namespace Rcpp;
// genEUTcpp
arma::mat genEUTcpp(arma::vec par, List dempars, int N, arma::mat A, arma::mat B, arma::mat pA, arma::mat pB, arma::vec Min, arma::vec Max, int ufunc, int lo, int unif, int choice_dbug);
RcppExport SEXP rcgen_genEUTcpp(SEXP parSEXP, SEXP demparsSEXP, SEXP NSEXP, SEXP ASEXP, SEXP BSEXP, SEXP pASEXP, SEXP pBSEXP, SEXP MinSEXP, SEXP MaxSEXP, SEXP ufuncSEXP, SEXP loSEXP, SEXP unifSEXP, SEXP choice_dbugSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< arma::vec >::type par(parSEXP);
Rcpp::traits::input_parameter< List >::type dempars(demparsSEXP);
Rcpp::traits::input_parameter< int >::type N(NSEXP);
Rcpp::traits::input_parameter< arma::mat >::type A(ASEXP);
Rcpp::traits::input_parameter< arma::mat >::type B(BSEXP);
Rcpp::traits::input_parameter< arma::mat >::type pA(pASEXP);
Rcpp::traits::input_parameter< arma::mat >::type pB(pBSEXP);
Rcpp::traits::input_parameter< arma::vec >::type Min(MinSEXP);
Rcpp::traits::input_parameter< arma::vec >::type Max(MaxSEXP);
Rcpp::traits::input_parameter< int >::type ufunc(ufuncSEXP);
Rcpp::traits::input_parameter< int >::type lo(loSEXP);
Rcpp::traits::input_parameter< int >::type unif(unifSEXP);
Rcpp::traits::input_parameter< int >::type choice_dbug(choice_dbugSEXP);
rcpp_result_gen = Rcpp::wrap(genEUTcpp(par, dempars, N, A, B, pA, pB, Min, Max, ufunc, lo, unif, choice_dbug));
return rcpp_result_gen;
END_RCPP
}
// genRDU
arma::mat genRDU(arma::vec par, List dempars, int N, arma::mat A, arma::mat B, arma::mat pA, arma::mat pB, arma::vec Min, arma::vec Max, int type, int ufunc, int unif, int lo, int choice_dbug);
RcppExport SEXP rcgen_genRDU(SEXP parSEXP, SEXP demparsSEXP, SEXP NSEXP, SEXP ASEXP, SEXP BSEXP, SEXP pASEXP, SEXP pBSEXP, SEXP MinSEXP, SEXP MaxSEXP, SEXP typeSEXP, SEXP ufuncSEXP, SEXP unifSEXP, SEXP loSEXP, SEXP choice_dbugSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< arma::vec >::type par(parSEXP);
Rcpp::traits::input_parameter< List >::type dempars(demparsSEXP);
Rcpp::traits::input_parameter< int >::type N(NSEXP);
Rcpp::traits::input_parameter< arma::mat >::type A(ASEXP);
Rcpp::traits::input_parameter< arma::mat >::type B(BSEXP);
Rcpp::traits::input_parameter< arma::mat >::type pA(pASEXP);
Rcpp::traits::input_parameter< arma::mat >::type pB(pBSEXP);
Rcpp::traits::input_parameter< arma::vec >::type Min(MinSEXP);
Rcpp::traits::input_parameter< arma::vec >::type Max(MaxSEXP);
Rcpp::traits::input_parameter< int >::type type(typeSEXP);
Rcpp::traits::input_parameter< int >::type ufunc(ufuncSEXP);
Rcpp::traits::input_parameter< int >::type unif(unifSEXP);
Rcpp::traits::input_parameter< int >::type lo(loSEXP);
Rcpp::traits::input_parameter< int >::type choice_dbug(choice_dbugSEXP);
rcpp_result_gen = Rcpp::wrap(genRDU(par, dempars, N, A, B, pA, pB, Min, Max, type, ufunc, unif, lo, choice_dbug));
return rcpp_result_gen;
END_RCPP
}
static const R_CallMethodDef CallEntries[] = {
{"rcgen_genEUTcpp", (DL_FUNC) &rcgen_genEUTcpp, 13},
{"rcgen_genRDU", (DL_FUNC) &rcgen_genRDU, 14},
{NULL, NULL, 0}
};
RcppExport void R_init_rcgen(DllInfo *dll) {
R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
R_useDynamicSymbols(dll, FALSE);
}

147
src/genCommon.cpp Normal file
View File

@ -0,0 +1,147 @@
#include <RcppArmadillo.h>
#include "genCommon.hpp"
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
// Function to add covariates
NumericVector getPars(int csize, List covv, List covi, NumericVector fpar) {
int pnum = 0;
double ii;
NumericVector index, var;
NumericVector par(csize);
for (int i = 0 ; i < csize ; i++) {
par[i] = fpar[pnum];
if ( covi[i] != R_NilValue) {
index = covi[i];
for (int j=0; j < index.size(); j++) {
pnum++;
ii = index[j];
var = covv[ii];
par[i] = par[i] + fpar[pnum] * var[0];
}
}
pnum++;
}
return(par);
}
// Functions needed for correlation
arma::mat fillcorr(arma::vec rho) {
int z = .5* ( sqrt(8*rho.n_elem + 1) + 1);
arma::mat RHO(z,z,arma::fill::eye);
int count = -1;
for (int row = 0; row < (z-1); row++) {
for (int col = row+1; col < z; col++) {
count++;
RHO(row,col) = rho[count];
}
}
RHO = symmatu(RHO);
return(RHO);
}
// Vector Power
arma::mat vpow(const arma::vec base, const arma::vec exp) {
arma::vec out(base.size());
std::transform(base.begin(), base.end(),
exp.begin(), out.begin(), ::pow);
return out;
}
// Multivariate normal inversion.
// Pass in unfiform distributed matrix, vector of means, and covaiance matirx, get back
// correlated uniform distirbuted matrix
NumericMatrix mvrnormInv(int n, arma::vec mu, arma::mat sigma) {
int ncols = sigma.n_cols;
arma::mat Y = arma::randn(n, ncols);
arma::mat RES = arma::repmat(mu, 1, n).t() + Y * arma::chol(sigma);
NumericMatrix OUT(n, ncols);
arma::colvec vtemp;
NumericVector htemp;
double mean;
double sd;
for (int i = 0; i < ncols ; i++) {
vtemp = RES.col(i);
htemp = as<NumericVector>(wrap(vtemp));
mean = mu[i];
sd = sqrt(sigma(i,i));
OUT(_,i) = pnorm( htemp , mean , sd );
}
return(OUT);
}
// Multivariate normal.
// Pass in unfiform distributed matrix, vector of means, and covaiance matirx, get back
// correlated normal distirbuted matrix
arma::mat mvrnorm(int n, arma::vec mu, arma::mat sigma) {
int ncols = sigma.n_cols;
arma::mat Y = arma::randn(n, ncols);
arma::mat RES = arma::repmat(mu, 1, n).t() + Y * arma::chol(sigma);
return(RES);
}
// Utility functions
arma::vec crra(arma::vec x, arma::vec r) {
return(vpow(x, 1 - r) / 1- r);
}
arma::mat corPrefs(arma::vec means, arma::vec sd, arma::vec rho) {
// How big is correlation matrix - zxz
int z = means.n_elem;
arma::mat RHO = fillcorr(rho);
arma::mat sigma(z,z,arma::fill::eye);
// Fill out the covariance matrix
for (int row = 0; row < z ; row++) {
for (int col = 0; col < z ; col++) {
sigma(row,col) = sd(row) * sd(col) * RHO(row,col);
}
}
// Get correlated normally distributed numbers
NumericVector norms = rnorm(z);
arma::mat Y(1, z);
Y.row(0) = as<arma::rowvec>(norms);
arma::mat out = means.t() + Y * arma::chol(sigma);
return(out);
}
arma::mat getUnif(arma::vec maxes, arma::vec mins) {
// How big is correlation matrix - zxz
int z = mins.n_elem;
NumericVector min2 = as<NumericVector>(wrap(mins));
NumericVector max2 = as<NumericVector>(wrap(maxes));
arma::vec means(z);
means.fill(0);
arma::mat sigma(z,z,arma::fill::eye);
// Get correlated normally distributed numbers
NumericVector norms = rnorm(z);
arma::mat Y(1, z);
Y.row(0) = as<arma::rowvec>(norms);
arma::mat normaldist = means.t() + Y * arma::chol(sigma);
NumericMatrix out = as<NumericMatrix>(wrap(normaldist));
out.row(0) = pnorm(out.row(0));
out.row(0) = out.row(0) * (max2 - min2) + min2;
return(as<arma::mat>(out));
}
double LN_mean(double m, double s) {
double out = log(m / sqrt(1 + ( pow(s, 2) / pow(m, 2) )));
return(out);
}
double LN_sd(double m, double s) {
double out = sqrt( log(1 + ( pow(s, 2) / pow(m, 2) )));
return(out);
}

31
src/genCommon.hpp Normal file
View File

@ -0,0 +1,31 @@
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
// Function to add covariates
NumericVector getPars(int csize, List covv, List covi, NumericVector fpar);
// Functions needed for correlation
arma::mat fillcorr(arma::vec rho);
// Vector Power
arma::mat vpow(const arma::vec base, const arma::vec exp);
// Multivariate normal inversion.
// Pass in unfiform distributed matrix, vector of means, and covaiance matirx, get back
// correlated normal distirbuted matrix
NumericMatrix mvrnormInv(int n, arma::vec mu, arma::mat sigma);
arma::mat mvrnorm(int n, arma::vec mu, arma::mat sigma);
// Utility functions
arma::vec crra(arma::vec x, arma::vec r);
// Get correlated preferences
arma::mat corPrefs(arma::vec means, arma::vec sd, arma::vec rho);
arma::mat getUnif(arma::vec mins, arma::vec maxes);
// Functions to transform means and standard deviations from normal to log normal
double LN_mean(double m, double s);
double LN_sd(double m, double s);

268
src/genEUT.cpp Normal file
View File

@ -0,0 +1,268 @@
#include <RcppArmadillo.h>
#include "genCommon.hpp"
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
// [[Rcpp::export]]
arma::mat genEUTcpp(arma::vec par, List dempars, int N, arma::mat A, arma::mat B,
arma::mat pA, arma::mat pB, arma::vec Min, arma::vec Max,
int ufunc, int lo, int unif, int choice_dbug) {
//int dbug=0;
//Rcout << "Here: " << ++dbug << std::endl;
// How big is the instrument
arma::uword isize = A.n_rows;
// How many columns of the final matrix are needed
int cnum = A.n_cols + pA.n_cols + B.n_cols + pB.n_cols + 4 + dempars.size();
// How many rows
int rnum = isize * N;
// How many outcome variables
int osize = A.n_cols;
// Per subject matix
double rm = par[0]; // CRRA Mean
double rs = par[1]; // CRRA Standard Deviation
double um = par[2]; // Fechner Mean
double us = par[3]; // Fechner Standard Deviation
double ru;
if (unif == 0) {
ru = par[4]; // Correlation CRRA and Fechner
}
else {
ru = 0; // Just a place holder
}
// columns containing individual distirbutional parameters
arma::mat rmmat(N,1);
arma::mat rsmat(N,1);
arma::mat ummat(N,1);
arma::mat usmat(N,1);
// Make the instrument matrix
arma::mat inst = join_rows(A, pA);
inst = join_rows(inst, B);
inst = join_rows(inst, pB);
inst = join_rows(inst, Max);
inst = join_rows(inst, Min);
// columns to be bound to the instrument
arma::mat rmat(isize,1);
arma::mat mumat(isize,1);
arma::mat idmat(isize,1);
arma::mat prefs(1,3);
double rmean = rm;
double rstd = rs;
if (lo == 1) {
rmean = LN_mean(1 - rm, rs);
rstd = LN_sd(1 - rm, us);
ru = -1 * ru;
}
double umean = LN_mean(um, us);
double ustd = LN_sd(um, us);
arma::vec means;
arma::vec sd;
arma::vec rho;
if (unif == 0) {
means = {rmean, umean};
sd = {rstd, ustd};
rho = {ru};
} else {
means = {rm, um};
sd = {rs, us};
rho = {ru};
}
// Declare utility vars
// Probability of A
arma::vec AProb(isize);
// For utilities
arma::mat uA = A;
arma::mat uB = A;
arma::mat UA(isize, 1);
arma::mat UB(isize, 1);
arma::mat UMin(isize, 1);
arma::mat UMax(isize, 1);
arma::mat CTX(isize, 1);
arma::mat PA(isize, 1);
arma::mat randomizer(isize, 1);
arma::mat choice(isize, 1);
// First subject
if (unif == 0) {
prefs = corPrefs(means, sd, rho);
} else {
prefs = getUnif(means, sd);
}
double r, mu;
r = prefs.at(0,0);
mu = prefs.at(0,1);
if (unif == 0) {
mu = exp(mu);
}
if (lo == 1) {
r = 1 - exp(r);
}
if (ufunc == 0) {
UMin = pow(Min, 1-r) / (1-r);
UMax = pow(Max, 1-r) / (1-r);
uA = pow(A, 1-r) / (1-r);
uB = pow(B, 1-r) / (1-r);
} else if (ufunc == 1) {
UMin = pow(Min, r);
UMax = pow(Max, r);
uA = pow(A, r);
uB = pow(B, r);
}
CTX = UMax - UMin;
// weight probabilities
uA %= pA;
uB %= pB;
UA = sum(uA,1);
UB = sum(uB,1);
// Rebase so UA = 0, add in context and fechner error
UB = (UB-UA) / CTX / mu;
// Probability of A
PA = 1.0 / (1.0 + exp(UB));
randomizer = arma::randu<arma::vec>(isize);
// If the probabiliy of A is greater than the number drawn from a uniform
// distribution, the choice is A otherwise it is B
std::transform(PA.begin(), PA.end(), randomizer.begin(), choice.begin(),
[] (double prob, double rand) {
if (prob > rand) return(0);
else return(1);
});
rmat.fill(r);
mumat.fill(mu);
idmat.fill(1);
// Attach ID
arma::mat per = join_rows(idmat, inst);
// If debugging, attach everything
if (choice_dbug == 1) {
per = join_rows( per, UA );
per = join_rows( per, UB );
per = join_rows( per, UMin );
per = join_rows( per, UMax );
per = join_rows( per, CTX );
per = join_rows( per, PA );
per = join_rows( per, 1 - PA );
per = join_rows( per, randomizer );
}
// Attach Choices
per = join_rows(per, choice);
// Attach preferences
per = join_rows(per, rmat);
per = join_rows(per, mumat);
// Output matirx
arma::mat out = per;
for (int i = 1; i < N; i++) {
// First subject
if (unif == 0) {
prefs = corPrefs(means, sd, rho);
} else {
prefs = getUnif(means, sd);
}
r = prefs.at(0,0);
mu = prefs.at(0,1);
if (lo == 1) {
r = 1 - exp(r);
}
if (unif == 0) {
mu = exp(mu);
}
if (ufunc == 0) {
UMin = pow(Min, 1-r) / (1-r);
UMax = pow(Max, 1-r) / (1-r);
uA = pow(A, 1-r) / (1-r);
uB = pow(B, 1-r) / (1-r);
} else if (ufunc == 1) {
UMin = pow(Min, r);
UMax = pow(Max, r);
uA = pow(A, r);
uB = pow(B, r);
}
CTX = UMax - UMin;
uA %= pA;
uB %= pB;
UA = sum(uA,1);
UB = sum(uB,1);
// Rebase so UA = 0, add in context and fechner error
UB = (UB-UA) / CTX / mu;
// Probability of A
PA = 1.0 / (1.0 + exp(UB));
// A random variable
randomizer = arma::randu<arma::vec>(isize);
// If the probabiliy of A is greater than the number drawn from a uniform
// distribution, the choice is A otherwise it is B
std::transform(PA.begin(), PA.end(), randomizer.begin(), choice.begin(),
[](double prob, double rand){
if (prob > rand) return(0);
else return(1);
});
// Attach stuff
rmat.fill(r);
mumat.fill(mu);
idmat.fill(i+1);
// Attach ID and instrument to each other
per = join_rows(idmat, inst);
// If debugging, attach everything
if (choice_dbug == 1) {
per = join_rows( per, UA );
per = join_rows( per, UB );
per = join_rows( per, UMin );
per = join_rows( per, UMax );
per = join_rows( per, CTX );
per = join_rows( per, PA );
per = join_rows( per, 1 - PA );
per = join_rows( per, randomizer );
}
// Attach Choices
per = join_rows(per, choice);
// Attach preferences
per = join_rows(per, rmat);
per = join_rows(per, mumat);
// Attach the to output matrix
out = join_cols(out, per);
}
return(out);
}

431
src/genRDU.cpp Normal file
View File

@ -0,0 +1,431 @@
#include <RcppArmadillo.h>
#include "genCommon.hpp"
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
// [[Rcpp::export]]
arma::mat genRDU(arma::vec par, List dempars, int N, arma::mat A, arma::mat B,
arma::mat pA, arma::mat pB, arma::vec Min, arma::vec Max,
int type, int ufunc, int unif, int lo, int choice_dbug) {
//int dbug=0;
// How big is the instrument
arma::uword isize = A.n_rows;
// How many columns of the final matrix are needed
int cnum = A.n_cols + pA.n_cols + B.n_cols + pB.n_cols + 4 + dempars.size();
// How many rows
int rnum = isize * N;
// How many outcome variables
int osize = A.n_cols;
arma::mat out(rnum,cnum);
double rm = par[0]; // CRRA Mean
double rs = par[1]; // CRRA Standard Deviation
double am = par[2]; // Alpha Mean
double as = par[3]; // Alpha Standard Deviation
// declare these things regardless of type
double bm, bs, um, us, ra, rb, ab, ru, au, bu;
if (type == 2) {
bm = par[4]; // Beta Mean
bs = par[5]; // Beta Standard Deviation
um = par[6]; // Fechner Mean
us = par[7]; // Fechner Standard Deviation
ra = par[8]; // Correlation CRRA and Alpha
rb = par[9]; // Correlation CRRA and Beta
ab = par[10]; // Correlation Alpha and Beta
ru = par[11]; // Correlation CRRA and Fechner
au = par[12]; // Correlation Alpha and Fechner
bu = par[13]; // Correlation Beta and Fechner
}
else {
um = par[4]; // Fechner Mean
us = par[5]; // Fechner Standard Deviation
ru = par[6]; // Correlation CRRA and Fechner
ra = par[7]; // Correlation CRRA and Alpha
au = par[8]; // Correlation Fechner and Alpha
}
if (unif == 1) {
ra = 0; // Correlation CRRA and Alpha
rb = 0; // Correlation CRRA and Beta
ab = 0; // Correlation Alpha and Beta
ru = 0; // Correlation CRRA and Fechner
au = 0; // Correlation Alpha and Fechner
bu = 0; // Correlation Beta and Fechner
}
// Make the instrument matrix
arma::mat inst = join_rows(A, pA);
inst = join_rows(inst, B);
inst = join_rows(inst, pB);
inst = join_rows(inst, Max);
inst = join_rows(inst, Min);
// columns to be bound to the instrument
arma::mat rmat(isize,1);
arma::mat amat(isize,1);
arma::mat bmat(isize,1);
arma::mat mumat(isize,1);
arma::mat idmat(isize,1);
arma::mat prefs(1,4);
// Sometimes we bound r to be less than 1
double rmean = rm;
double rstd = rs;
if (lo == 1) {
rmean = LN_mean(1 - rm, rs);
rstd = LN_sd(1 - rm, us);
ru = -1 * ru;
ra = -1 * ra;
if (type ==2) {
rb = -1 * rb;
}
}
// some parameters need to have their distributions bound to >0
double umean = LN_mean(um, us);
double ustd = LN_sd(um, us);
double amean = LN_mean(am, as);
double astd = LN_sd(am, as);
arma::vec means, sd, rho;
if (type ==2) {
double bmean = LN_mean(bm, bs);
double bstd = LN_sd(bm, bs);
means = {rmean, amean, bmean, umean};
sd = {rstd, astd, bstd, ustd};
rho = {ra, rb, ab, ru, au, bu};
if (unif == 1) {
means = {rm, am, bm, um};
sd = {rs, as, bs, us};
rho = {ra, rb, ab, ru, au, bu};
}
} else {
means = {rm, amean, umean};
sd = {rs, astd, ustd};
rho = {ra, ru, au};
if (unif == 1) {
means = {rm, am, um};
sd = {rs, as, us};
rho = {ra, ru, au};
}
}
// Declare utility vars
// Probability of A
arma::vec AProb(isize);
// For utilities
arma::mat uA = A;
arma::mat uB = A;
arma::vec UA(isize);
arma::vec UB(isize);
arma::vec UMin(isize);
arma::vec UMax(isize);
arma::vec CTX(isize);
arma::vec PA(isize);
arma::vec randomizer(isize);
arma::mat choice(isize,1);
// Cummulative probabilities
arma::mat pA_f = arma::cumsum(pA,1);
// Make a matrix that doesn't have the the lowest probabiliy choice
arma::mat pAc_f = pA_f.cols(0, pA_f.n_cols - 2);
// Add a column at the begining, filled with 0s
pAc_f.insert_cols(0,1);
arma::mat pB_f = arma::cumsum(pB,1);
arma::mat pBc_f = pB_f.cols(0, pB_f.n_cols - 2);
pBc_f.insert_cols(0,1);
// For decision weights
arma::mat dA = pA;
arma::mat dB = pB;
arma::mat dAc = pA;
arma::mat dBc = pB;
// First subject
if (unif == 0) {
prefs = corPrefs(means, sd, rho);
} else {
prefs = getUnif(means, sd);
}
double r, alpha, beta, mu;
r = prefs.at(0,0);
if (lo == 1) {
r = 1 - exp(r);
}
alpha = prefs.at(0,1);
if (type == 2) {
beta = prefs.at(0,2);
mu = prefs.at(0,3);
} else {
mu = prefs.at(0,2);
}
if (unif == 0) {
alpha = exp(alpha);
mu = exp(mu);
if (type == 2) {
beta = exp(beta);
}
}
if (ufunc == 0) {
UMin = pow(Min, 1.0-r) / (1.0-r);
UMax = pow(Max, 1.0-r) / (1.0-r);
uA = pow(A, 1.0-r) / (1.0-r);
uB = pow(B, 1.0-r) / (1.0-r);
} else if (ufunc == 1) {
UMin = pow(Min, r);
UMax = pow(Max, r);
uA = pow(A, r);
uB = pow(B, r);
}
CTX = UMax - UMin;
// Add probability weighting - this function uses power function
if (type == 0) {
dA = pow(pA_f, alpha);
dAc = pow(pAc_f, alpha);
dA -= dAc;
dB = pow(pB_f, alpha);
dBc = pow(pBc_f, alpha);
dB -= dBc;
} else if (type == 1) {
dA = pow(pA_f, alpha) / pow(pow(pA_f, alpha) + pow(1-pA_f, alpha), 1/alpha);
dAc = pow(pAc_f, alpha) / pow(pow(pAc_f, alpha) + pow(1-pAc_f, alpha), 1/alpha);
dA -= dAc;
dB = pow(pB_f, alpha) / pow(pow(pB_f, alpha) + pow(1-pB_f, alpha), 1/alpha);
dBc = pow(pBc_f, alpha) / pow(pow(pBc_f, alpha) + pow(1-pBc_f, alpha), 1/alpha);
dB -= dBc;
} else if (type == 2) {
dA = exp(-beta * (pow(-log(pA_f), alpha)));
dAc = exp(-beta * (pow(-log(pAc_f), alpha)));
dA.elem( find_nonfinite(dA) ).zeros();
dAc.elem( find_nonfinite(dAc) ).zeros();
dA -= dAc;
dA.elem( find_nonfinite(dA) ).zeros();
dB = exp(-beta * (pow(-log(pB_f), alpha)));
dBc = exp(-beta * (pow(-log(pBc_f), alpha)));
dB.elem( find_nonfinite(dB) ).zeros();
dBc.elem( find_nonfinite(dBc) ).zeros();
dB -= dBc;
dB.elem( find_nonfinite(dB) ).zeros();
}
uA %= dA;
uB %= dB;
UA = sum(uA,1);
UB = sum(uB,1);
// Rebase so UA = 0, add in context and fechner error
UB = (UB-UA) / CTX / mu;
// Probability of A
PA = 1.0 / (1.0 + exp(UB));
randomizer = arma::randu<arma::vec>(isize);
// If the probabiliy of A is greater than the number drawn from a uniform
// distribution, the choice is A otherwise it is B
std::transform(PA.begin(), PA.end(), randomizer.begin(), choice.begin(),
[](double prob, double rand){
if (prob > rand) return(0);
else return(1);
});
rmat.fill(r);
amat.fill(alpha);
if (type == 2) bmat.fill(beta);
mumat.fill(mu);
idmat.fill(1);
// Attach ID
arma::mat per = join_rows(idmat, inst);
// If debugging, attach everything
if (choice_dbug == 1) {
per = join_rows( per, dA );
per = join_rows( per, dB );
per = join_rows( per, UA );
per = join_rows( per, UB );
per = join_rows( per, UMin );
per = join_rows( per, UMax );
per = join_rows( per, CTX );
per = join_rows( per, PA );
per = join_rows( per, 1 - PA );
per = join_rows( per, randomizer );
}
// Attach Choices
per = join_rows(per, choice);
// Attach preferences
per = join_rows(per, rmat);
per = join_rows(per, amat);
if (type == 2) per = join_rows(per, bmat);
per = join_rows(per, mumat);
// Output matirx
out = per;
for (int i =1; i < N; i++) {
// First subject
if (unif == 0) {
prefs = corPrefs(means, sd, rho);
} else {
prefs = getUnif(means, sd);
}
r = prefs.at(0,0);
if (lo == 1) {
r = 1 - exp(r);
}
alpha = prefs.at(0,1);
if (type == 2) {
beta = prefs.at(0,2);
mu = prefs.at(0,3);
} else {
mu = prefs.at(0,2);
}
if (unif == 0) {
alpha = exp(alpha);
mu = exp(mu);
if (type == 2) {
beta = exp(beta);
}
}
if (ufunc == 0) {
UMin = pow(Min, 1-r) / (1-r);
UMax = pow(Max, 1-r) / (1-r);
uA = pow(A, 1-r) / (1-r);
uB = pow(B, 1-r) / (1-r);
} else if (ufunc == 1) {
UMin = pow(Min, r);
UMax = pow(Max, r);
uA = pow(A, r);
uB = pow(B, r);
}
CTX = UMax - UMin;
// Add probability weighting - 0 for power, 1 for inverse-S
if (type == 0) {
dA = pow(pA_f, alpha);
dAc = pow(pAc_f, alpha);
dA -= dAc;
dB = pow(pB_f, alpha);
dBc = pow(pBc_f, alpha);
dB -= dBc;
} else if (type == 1) {
dA = pow(pA_f, alpha) / pow(pow(pA_f, alpha) + pow(1-pA_f, alpha), 1/alpha);
dAc = pow(pAc_f, alpha) / pow(pow(pAc_f, alpha) + pow(1-pAc_f, alpha), 1/alpha);
dA -= dAc;
dB = pow(pB_f, alpha) / pow(pow(pB_f, alpha) + pow(1-pB_f, alpha), 1/alpha);
dBc = pow(pBc_f, alpha) / pow(pow(pBc_f, alpha) + pow(1-pBc_f, alpha), 1/alpha);
dB -= dBc;
} else if (type == 2) {
dA = exp(-beta * (pow(-log(pA_f), alpha)));
dAc = exp(-beta * (pow(-log(pAc_f), alpha)));
dA.elem( find_nonfinite(dA) ).zeros();
dAc.elem( find_nonfinite(dAc) ).zeros();
dA -= dAc;
dA.elem( find_nonfinite(dA) ).zeros();
dB = exp(-beta * (pow(-log(pB_f), alpha)));
dBc = exp(-beta * (pow(-log(pBc_f), alpha)));
dB.elem( find_nonfinite(dB) ).zeros();
dBc.elem( find_nonfinite(dBc) ).zeros();
dB -= dBc;
dB.elem( find_nonfinite(dB) ).zeros();
}
uA %= dA;
uB %= dB;
UA = sum(uA,1);
UB = sum(uB,1);
// Rebase so UA = 0, add in context and fechner error
UB = (UB-UA) / CTX / mu;
// Probability of A
PA = 1.0 / (1.0 + exp(UB));
// A random variable
randomizer = arma::randu<arma::vec>(isize);
// If the probabiliy of A is greater than the number drawn from a uniform
// distribution, the choice is A otherwise it is B
std::transform(PA.begin(), PA.end(), randomizer.begin(), choice.begin(),
[] (double prob, double rand) {
if (prob > rand) return(0.0);
else return(1.0);
});
// Attach stuff
rmat.fill(r);
amat.fill(alpha);
if (type == 2) bmat.fill(beta);
mumat.fill(mu);
idmat.fill(i + 1);
// Attach ID
arma::mat per = join_rows(idmat, inst);
// If debugging, attach everything
if (choice_dbug == 1) {
per = join_rows( per, dA );
per = join_rows( per, dB );
per = join_rows( per, UA );
per = join_rows( per, UB );
per = join_rows( per, UMin );
per = join_rows( per, UMax );
per = join_rows( per, CTX );
per = join_rows( per, PA );
per = join_rows( per, 1 - PA );
per = join_rows( per, randomizer );
}
// Attach Choices
per = join_rows(per, choice);
// Attach preferences
per = join_rows(per, rmat);
per = join_rows(per, amat);
if (type == 2) per = join_rows(per, bmat);
per = join_rows(per, mumat);
// Attach the to output matrix
out = join_cols(out, per);
}
return(out);
}