##
## Copyright 2016 Brent Kaplan
##
## This file is part of beezdemand.
##
## beezdemand is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, version 2.
##
## beezdemand is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with beezdemand. If not, see <http://www.gnu.org/licenses/gpl-2.0.html>.
##
## summary
## R script for simulating data
##
##' Simulate demand data
##'
##' Generates and saves simulated datasets in the manner specified in Koffarnus, Franck, Stein, & Bickel (2015).
##' @title Simulate Demand Data
##' @param nruns Number of runs. Default value is 10
##' @param setparams A 6x1 matrix (or 6 element vector) containing (in order) mean log10alpha, sd log10alpha, mean log10q0, sd log10q0, k, sd of consumption values across all prices
##' @param sdindex A vector of n length of sd consumption values for n prices
##' @param x A vector of n prices
##' @param outdir Optional. Directory to save results. Must end with a "/"
##' @param fn Optional. Filename of saved RData object
##' @return Invisibly a list consisting of: rounded consumption values, unrounded consumption values, simulation parameters, and inState and outState of seeds.
##' @author Brent Kaplan <bkaplan.ku@@gmail.com>
##' @examples
##' ## set values
##' setparams <- vector(length = 4)
##' setparams <- c(-2.5547, .702521, 1.239893, .320221, 3.096, 1.438231)
##' names(setparams) <- c("alphalm", "alphalsd", "q0lm", "q0lsd", "k", "yvalssd")
##' sdindex <- c(2.1978, 1.9243, 1.5804, 1.2465, 0.8104, 0.1751, 0.0380, 0.0270)
##' x <- c(.1, 1, 3, 10, 30, 100, 300, 1000)
##' set.seed(1234)
##' sim <- SimulateDemand(nruns = 1, setparams = setparams, sdindex = sdindex, x = x)
##' sim
##' @export
SimulateDemand <- function(nruns = 10, setparams, sdindex, x, outdir = NULL, fn = NULL) {
RunOneSim <- function(run = 1, setparams, sdindex, x) {
## Save inSate seed
inState <- .Random.seed
## Initialize vector for simulation parameters
simparams <- vector(length = 4)
## Store simulation parameters and label
simparams[1] <- rnorm(1, setparams[1], setparams[2])
simparams[2] <- rnorm(1, setparams[3], setparams[4])
simparams[3] <- 10^simparams[1]
simparams[4] <- 10^simparams[2]
names(simparams) <- c("alphalr", "q0lr", "alphar", "q0r")
## Save outState seed
outState <- .Random.seed
## Calculate consumption values at each price
sim <- function(x, sdindex) {
m <- 10^(simparams["q0lr"] + setparams[5] *
((exp((-1 * simparams["alphar"]) * simparams["q0r"] * x)) - 1))
s <- sdindex * setparams[6]
rnorm(1, m, s)
}
ysim <- mapply(sim, x, sdindex)
## Round resulting consumption values, truncating at 0
ysimr <- pmax(round(ysim, 0), 0)
## Return list: unrounded and rounded consumption values,
## simulation parameters, inState and outState of seed
list(dat = data.frame("y" = ysim, "yr" = ysimr),
"simparams" = simparams, "inState" = inState, "outState" = outState)
}
runs <- seq_len(nruns)
manysims <- lapply(runs, RunOneSim, setparams, sdindex, x)
simvalues <- sapply(manysims, "[[", 1)
simparams <- sapply(manysims, "[[", 2)
seeds <- sapply(manysims, "[", 3:4)
simvaluesraw <- simvalues[-2, , drop = FALSE]
simvalues <- simvalues[-1, , drop = FALSE]
simvaluesraw <- do.call("cbind", simvaluesraw)
simvalues <- do.call("cbind", simvalues)
simvaluesraw <- cbind(x, simvaluesraw)
simvalues <- cbind(x, simvalues)
colnames(simvaluesraw) <- c("x", runs)
colnames(simvalues) <- c("x", runs)
rownames(simparams) <- c("alphalr", "q0lr", "alphar", "q0r")
colnames(simparams) <- runs
simparams <- t(simparams)
colnames(seeds) <- runs
seeds <- t(seeds)
## Reshape long
simvalues <- data.frame(simvalues)
simvalues <- reshape2::melt(simvalues, id.vars = "x")
colnames(simvalues) <- c("x", "id", "y")
simvalues <- simvalues[, c("id", "x", "y")]
simvalues$id <- gsub("X", "", simvalues$id)
simvalues$id <- as.numeric(simvalues$id)
simvaluesraw <- data.frame(simvaluesraw)
simvaluesraw <- reshape2::melt(simvaluesraw, id.vars = "x")
colnames(simvaluesraw) <- c("x", "id", "y")
simvaluesraw <- simvaluesraw[, c("id", "x", "y")]
simvaluesraw$id <- gsub("X", "", simvaluesraw$id)
simvaluesraw$id <- as.numeric(simvaluesraw$id)
if (!is.null(outdir) && !is.null(fn)) {
if (!dir.exists(outdir)) {
dir.create(outdir)
}
saveobjects <- c("simvalues", "simvaluesraw", "simparams", "seeds")
save(list = saveobjects, file = paste0(outdir, fn, ".RData"))
}
invisible(list("simvalues" = simvalues, "simvaluesraw" = simvaluesraw,
"simparams" = simparams, "seeds" = seeds))
}
##' Gets values used in SimulateDemand
##'
##' Gets values used in SimulateDemand
##' @title Get Values for SimulateDemand
##' @param dat Dataframe (long form)
##' @return List of 3: setaparams, sdindex, x
##' @author Brent Kaplan <bkaplan.ku@@gmail.com>
##' @examples
##' GetValsForSim(apt)
##' @export
GetValsForSim <- function(dat) {
k <- log10(max(dat[dat$y > 0, "y"])) - log10(min(dat[dat$y > 0, "y"]))
## Get N unique participants, prices informing loops
participants <- unique(dat$id)
np <- length(participants)
prices <- unique(dat$x)
nprices <- length(prices)
cnames <- c("Participant", "Q0", "Alpha", as.numeric(prices))
dfres <- data.frame(matrix(vector(),
np,
length(cnames),
dimnames = list(c(), c(cnames))),
stringsAsFactors = FALSE)
for (i in seq_len(np)) {
dfres[i, "Participant"] <- participants[i]
adf <- NULL
adf <- dat[dat$id == participants[i], ]
adf[, "k"] <- k
fit <- NULL
fit <- try(nlsr::wrapnlsr(data = adf,
y ~ q0 * 10^(k * (exp(-alpha * q0 * x) - 1)),
start = list(q0 = 10, alpha = 0.01)),
silent = TRUE)
if (!inherits(fit, "try-error")) {
dfres[i, c("Q0", "Alpha")] <- as.numeric(coef(fit)[c("q0", "alpha")])
dfres[i, 4:NCOL(dfres)] <- resid(fit)
}
}
setparams <- vector("numeric", 6)
setparams[1] <- mean(log10(dfres$Alpha), na.rm = TRUE)
setparams[2] <- sd(log10(dfres$Alpha), na.rm = TRUE)
setparams[3] <- mean(log10(dfres$Q0), na.rm = TRUE)
setparams[4] <- sd(log10(dfres$Q0), na.rm = TRUE)
setparams[5] <- k
setparams[6] <- sd(dat$y)
names(setparams) <- c("alphalm", "alphalsd", "q0lm", "q0lsd", "k", "yvalssd")
sdindex <- vector("numeric", nprices)
for (i in 4:NCOL(dfres)) {
j <- i - 3
sdindex[j] <- sd(dfres[ , i], na.rm = TRUE)
}
list("setparams" = setparams, "sdindex" = sdindex, "x" = prices)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.