ericnovik

library(ggplot2)
library(deSolve)
library(tidyr)
theme_set(theme_minimal())

# I: infected, S: susceptible, R: recovered
sir <- function(t, state, pars) {
  with(as.list(c(state, pars)), {
    N <- S + I + R
    dSdt <- -b * I/N * S

ericnovik

library(dplyr)
library(tidyr)
library(magrittr)
library(bayesplot)

plot_areas <- function(d, n, mu, s, cutoff, left = TRUE, h, v) {
  
  if (left) {
    Pr <- apply(d, 2, function (x) mean(x < cutoff))
  } else {

ericnovik

library(manipulate)
manipulate({
  lambda <- seq(0, 5, length = 100)
  prior <- dgamma(lambda, a, b)
  post <- dgamma(lambda, a + sum(y), b + length(y))
  tibble(lambda, prior, post) %>%
    ggplot(aes(lambda, prior)) + geom_line() +
    geom_line(aes(y = post), color = "red") +
    ylab("density") +
    xlab(expression(lambda))

ericnovik

data {
  int<lower=2> K; // number of categories
  int<lower=0> N; // number of observations
  int<lower=1> D; // number of predictors
  int<lower=1, upper=K> y[N];
  matrix[N, D] x;
}
parameters {
  vector[D] beta;
  ordered[K-1] c; // cutpoints

ericnovik

library(cranlogs)
library(ggplot2)
d <- cran_downloads(
  package = c("rstan", "brms", "rstanarm"),
  from = "2015-01-01",
  to = Sys.Date()
)

ggplot(d, aes(x = date, y = count, color = package)) +
  geom_smooth() +

ericnovik

library(cowplot)
f <- function(x) {
  exp(-x^2) + sin(10*x)
}
n <- 3e3
bounds <- c(0, 1)
set.seed(123)
x <- runif(n, bounds[1], bounds[2])
y <- f(x)
p1 <- qplot(x, y, geom = "line")

ericnovik

library(cowplot)
library(tidyverse)

pow <- function(x, p) sign(x) * abs(x)^p
n <- 1e3; num_of_curves <- 10;
x <- seq(-3, 3, length.out = n)
C <- seq(-4, 4, length.out = num_of_curves)
C <- c(C, 0) # include the C = 0 solution
ncol <- length(C)
y <- matrix(nrow = n, ncol =ncol)

ericnovik

library(cowplot)
n <- 1e4
X <- rlogis(n)
Y <- plogis(X)
plot_dens <- function(data, ...) {
  qplot(
    data,
    geom = "histogram",
    alpha = I(1 / 2), 
    ...

ericnovik

HMC = function (U, grad_U, epsilon, L, current_q)
{
  q = current_q
  p = rnorm(length(q), 0, 1) # independent standard normal variates
  current_p = p
  # Make a half step for momentum at the beginning
  p = p - epsilon * grad_U(q) / 2
  # Alternate full steps for position and momentum
  for (i in 1:L)
  {

ericnovik

x <- matrix(rep(10, 12), nrow = 4)
y <- c(2, 5, 10)
div_matrix_by_vec <- function(mat, vec) {
  stopifnot(is.matrix(mat))
  stopifnot(is.numeric(vec))
  stopifnot(ncol(mat) == length(vec))
  
  t(apply(mat, 1, function(x) x/vec))
}
all.equal(div_matrix_by_vec(x, y),