Skip to content

Instantly share code, notes, and snippets.

View statwonk's full-sized avatar
🏠
Working from home

Christopher Peters statwonk

🏠
Working from home
View GitHub Profile
@statwonk
statwonk / prediction_intervals.R
Created June 25, 2019 22:49
An example of checking the statistical coverage properties of an algorithm.
library(tidyverse)
tibble(boot = seq_len(1e2)) %>%
mutate(prediction_intervals = map(boot, function(i) {
rpois(100, 10) -> x
(sum(x[6:100]) -> total_oracle)
# Given the first five, what can we say about the sum of 95 more?
x[1:5] -> first_five
@statwonk
statwonk / normals.R
Created June 19, 2019 04:27
Part of the Normal family of distributions.
library(tidyverse)
plot.new()
par(bg = "black")
plot(density(rnorm(1e6)), col = "deeppink", lwd = 1,
xlab = "", ylab = "", main = "Normal distribution", axes = FALSE)
seq(1, 10, 0.05) %>%
map(function(x){ rnorm(1e5, 0, x)}) %>%
map(function(x) { density(x) %>%
lines(col = colors() %>% Filter(function(.x) grepl("pink", .x), .) %>% sample(1),
lwd = sample(seq(0.1, 0.1, 0.01), 1),
@statwonk
statwonk / inezs_weibulls.R
Created June 1, 2019 22:06
;O{) # o-----\. # /.;/?}
library(tidyverse)
plot.new()
par(bg = "black")
plot(density(rweibull(1e6, 2, 0.1)), col = "pink3", lwd = 1, ylim = c(0, 20),
xlab = "", ylab = "", main = "Weibull distribution", axes = FALSE)
seq(1, 5, 0.25) %>%
map(function(x){ rweibull(1e6, x, 0.1)}) %>%
map(function(x) { density(x) %>%
lines(col = sample(colors(), 1),
lwd = sample(seq(0.1, 3, 0.1), 1),
@statwonk
statwonk / control-vs-interaction.R
Created May 17, 2019 17:12
A comparison of std.errors of the treatment effect from an additive DGP experiment where a pre-treatment variable is interacted maximally, only additive no interaction, or unadjusted for.
library(tidyverse)
library(lme4)
library(broom)
N <- 1e4
generate_data <- function() {
tibble(
intercept = 1,
pretreatment_var = case_when(rbinom(N, 1, p = 0.5) == 1 ~ 0.5, TRUE ~ -0.5),
@statwonk
statwonk / cdfs.R
Created April 11, 2019 18:08
CDF gist
library(tidyverse)
# Cumulative density functions
?hist
# built in datasets
str(mtcars)
head(mtcars)
@statwonk
statwonk / km-quantiles.R
Created April 7, 2019 13:21
A routine showing how to calculate quantiles.
library(survival)
library(Hmisc)
# 1. a weibull random number generator, see http://statwonk.com/weibull.html
rweibull_cens <- function(n, shape, scale) {
# will happen to see the death time first or censoring?
rweibull(n, shape = shape, scale = scale) -> a_random_death_time
rweibull(n, shape = shape, scale = scale) -> a_random_censor_time
pmin(a_random_censor_time, a_random_death_time) -> observed_time

Keybase proof

I hereby claim:

  • I am statwonk on github.
  • I am statwonk (https://keybase.io/statwonk) on keybase.
  • I have a public key ASBKY0yAKx1KghVEq7DmzhrFWAIZn94O-vhQHcjU9Wwz-go

To claim this, I am signing this object:

@statwonk
statwonk / exp_relative_likelihoods.R
Created October 23, 2017 04:14
100 exponential random variables and their relative likelihoods
suppressPackageStartupMessages({
library(dplyr)
library(ggplot2)
library(purrr)
})
exp_likelihood <- function(x, lambda) {
(lambda^length(x))*exp(-lambda*length(x)*mean(x))
}
seq_len(1e2) %>%
lapply(function(i) {
@statwonk
statwonk / library.R
Created July 6, 2016 02:31
R's library() function
function (package, help, pos = 2, lib.loc = NULL, character.only = FALSE,
logical.return = FALSE, warn.conflicts = TRUE, quietly = FALSE,
verbose = getOption("verbose"))
{
testRversion <- function(pkgInfo, pkgname, pkgpath) {
if (is.null(built <- pkgInfo$Built))
stop(gettextf("package %s has not been installed properly\n",
sQuote(pkgname)), call. = FALSE, domain = NA)
R_version_built_under <- as.numeric_version(built$R)
if (R_version_built_under < "3.0.0")
@statwonk
statwonk / .getRequiredPackages2.R
Created July 6, 2016 02:09
.getRequiredPackages2
function (pkgInfo, quietly = FALSE, lib.loc = NULL, useImports = FALSE)
{
.findVersion <- function(pkg, lib.loc = NULL) {
pfile <- system.file("Meta", "package.rds", package = pkg,
lib.loc = lib.loc)
if (nzchar(pfile))
as.numeric_version(readRDS(pfile)$DESCRIPTION["Version"])
else NULL
}
.findAllVersions <- function(pkg, lib.loc = NULL) {