jkeirstead

# Generate a Monte Carlo sample using Sobol' low-discrepancy quasi-random sequences
# James Keirstead
# 3 February 2012
#
# Random sampling with R's standard methods is inefficient for Monte Carlo analysis as
# the sampled values do not cover the parameter space evenly.  This Gist allows users
# to create parameter samples using Sobol' sequences to get around this problem.

# makeMCSample
#    Makes a Monte Carlo sample using Sobol' sequences

jkeirstead

# Demo of Gaussian process regression with R
# James Keirstead
# 5 April 2012

# Chapter 2 of Rasmussen and Williams's book `Gaussian Processes
# for Machine Learning' provides a detailed explanation of the
# math for Gaussian process regression.  It doesn't provide
# much in the way of code though.  This Gist is a brief demo
# of the basic elements of Gaussian process regression, as
# described on pages 13 to 16.

jkeirstead

## Analysis of DECC NEED data
## James Keirstead
## 22 November 2012
##
## Underlying data available from
##     http://www.decc.gov.uk/en/content/cms/statistics/energy_stats/en_effic_stats/need/need.aspx
## For more info, see
##     http://www.jameskeirstead.ca/blog/we-need-more-data/ 

jkeirstead

#' Calculate analytic hierarchy process weights
#'
#' The analytic hierarchy process (AHP) can be used for multi-criteria
#' decision analysis.  It performs a pairwise comparison of a list of
#' options, with a respect to a particular goal.
#'
#' @detail The script will ask you to compare pairs of options.
#'  Answer each question with respect to the statement goal and use
#'  scores from 1,3,5,7,9 (or their inverse) where 1 = equally
#'  important, 3 = moderately preferred, 5 = strongly preferred, 7 =

jkeirstead

## District heating AHP example
##
## This example demonstrates how the analytic hierarchy process can be
## used to select a preferred technique for forecasting the thermal
## load of a district heating system.
##
## James Keirstead (j.keirstead@imperial.ac.uk)
## Presented 25 June 2013 at ISIE Conference, Ulsan, South Korea
##

jkeirstead

// Predict from Gaussian Process
// All data parameters must be passed as a list to the Stan call
// Based on original file from https://code.google.com/p/stan/source/browse/src/models/misc/gaussian-process/

data {
int<lower=1> N1;
vector[N1] x1;
vector[N1] y1;
int<lower=1> N2;
vector[N2] x2;

jkeirstead

## List the functions in a file, or if not specified the global environment
##
## @param src an optional source file
##
## @return a character vector of function names
ls_functions <- function(src=NULL) {
  if (is.null(src)) {
    classes <-   sapply(ls(.GlobalEnv), function(x) class(eval(parse(text=x))))
  } else {
    e <- new.env()

jkeirstead

# Analysis of EPSRC fellowship success rates
# 18 November 2013
# James Keirstead

##' Calculates the bounds for a funnel plot
##'
##' Calculates the upper and lower confidence interval limits for use
##' in a funnel plot.  Assumes a binomial discrete distribution with a
##' probability of success theta.
##' 

jkeirstead

##' Converts a vector of values to a formatted prose list
##'
##' @param vals a vector of values
##' @param oxford a boolean indicating whether Oxford comma should be used
##' @param and word before ultimate entry
##' @return a character string
prose_vector <- function(vals, oxford=FALSE, and="and") {
    if (length(vals)>1) {
        start <- head(vals, -1)
        start <- paste0(start, collapse=", ")

jkeirstead

## Define the problem data
x <- read.csv("x-values.csv")$x

## Fit the power law. Looks okay with plots
library(poweRlaw)
mod <- conpl$new(x)    
xmin <- estimate_xmin(mod, xmins=head(sort(x), n=10))
mod$setXmin(xmin)
pars <- estimate_pars(mod)
mod$setPars(pars)