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 / beta_interval_data.R
Last active December 31, 2020 18:26
Using beta-distributed interval-censored data to produce an estimate for the median share of US adults having taken the vaccine by end of Q2 2021.
library(tidyverse)
library(fitdistrplus)
dplyr::select -> select
.Machine$double.eps -> eps
1975 -> N # number of responses
tibble(lower = c(0, 0.25, 0.5, 0.75), # lower bins
upper = c(0.25 + eps, 0.5 + eps, 0.75 + eps, 1), # upper bins
pct = c(0.32, 0.51, 0.15, 1 - sum(0.32, 0.51, 0.15)), # response shares
n = floor(pct * N)) %>% # implied responses + eps
@statwonk
statwonk / mktcap.R
Last active December 25, 2020 21:53
Some rough beliefs about market capitalization for a new entrant from a variety of sectors, sub-sectors and states.
library(tidyverse)
library(rvest)
library(gamlss)
library(brms)
library(tidybayes)
select <- dplyr::select
####################################################################################
# Model the market capitalizations of members of the S&P 500.
####################################################################################
@statwonk
statwonk / gist:6283c5b01e5896b94c46edfdd9ff490a
Created December 25, 2020 21:12
A model of the market capitalizations of S&P 500 members.
library(tidyverse)
library(rvest)
library(gamlss)
library(brms)
library(tidybayes)
select <- dplyr::select
####################################################################################
# Model the market capitalizations of members of the S&P 500.
####################################################################################
library(tidyverse)
library(quantmod)
library(gamlss)
select <- dplyr::select
posix <- function(x) { as.POSIXct(x, origin = "1970-01-01") }
## "Fat tails"
## Here we compare the residuals from the normal and t distributional models.
## Notice standardized error is worse in the normal model. This happens
## because returns are leptokurtic (large surprises should be expected, there's risk in stock returns),
@statwonk
statwonk / risk_adds_up2.R
Last active September 2, 2020 13:10
In the context of coronavirus, let's revisit this post about risk accumulation. https://twitter.com/statwonk/status/1160542394544267265?s=20
library(tidyverse)
library(ggthemes)
expand.grid(
risk = seq(0.1/5e3, 1/5e3, 1e-05), # average daily risk e.g. - 1,000 infected per day in Alabama / 5,000,000 AL population
units_of_exposure = seq_len(31) # days of exposure (up to 31 days)
) %>% as_tibble() %>%
mutate(total_risk = map2_dbl(risk, units_of_exposure, ~ 1 - (1 - .x)^(.y)),
total_odds = 1/total_risk,
risk_threshold = case_when(total_odds <= 5e2 ~ "Worse than 1 in 500",
total_odds <= 1e3 ~ "Worse than 1 in 1k chance",
@statwonk
statwonk / generate_multilevel_logistic_data.R
Created August 30, 2020 19:25
A tool to simulate multilevel logistic data.
library(tidyverse)
library(brms)
library(tidybayes)
1e7 -> N # obs
1 -> J # groups of members
10 -> K # members
0.5 -> base_p # base rate, this is logistic regression
# sample member coefficients
@statwonk
statwonk / out_of_core.py
Created August 30, 2020 14:24
scikit-learn partial_fit out-of-core learning.
import sklearn
from sklearn import naive_bayes
import pandas as pd
import numpy as np
d = pd.read_csv("data.csv")
y = d.iloc[:, 1]
X = d.iloc[:,list(range(2, d.shape[1]))]
@statwonk
statwonk / variational_logistic.R
Created August 15, 2020 21:54
A simulation to learn about variational inference and compare it to MCMC.
library(tidyverse)
library(brms)
library(tidybayes)
3e4 -> N
40 -> K
rnorm(K) -> group_coefs
tibble(K = factor(rep(paste0("group_", seq_len(K)), length.out = N))) %>%
mutate(coef = rep(group_coefs, N/40)) %>%
@statwonk
statwonk / coronavirus.R
Last active March 1, 2020 17:19
Showing that the number of deaths by coronavirus is somewhat robust to CFR restricted above 1% (South Korea lower bound at 0.4%).
set.seed(1)
N <- 1e7 # sims
quantiles_of_interest <- function(x) { quantile(x, c(0.0001, 0.25, 0.5, 0.75, 0.9999)) }
death_rate <- function() { pmin(pmax(rlnorm(N, log(0.02), 0.35), 0.001), 0.06) }
quantiles_of_interest(death_rate())
# very robust / ignorant belief of range of case 2.4k to 60M cases.
susceptible_cases <- function() { runif(N, 0.0024, 150) }
@statwonk
statwonk / coronavirus.R
Created March 1, 2020 16:29
Coronavirus estimates
set.seed(0)
N <- 1e7 # sims
quantiles_of_interest <- function(x) { quantile(x, c(0.0001, 0.25, 0.5, 0.75, 0.9999)) }
death_rate <- function() { pmin(pmax(rlnorm(N, log(0.02), 0.35), 0.001), 0.06) }
# Statwonk's belief of final US death rate
quantiles_of_interest(death_rate())
# 0.01% 25% 50% 75% 99.99%