lzachmann

# Moment match for the beta distribution.
params_from_moments <- function(mu, sigma) {
  alpha <- (mu^2 - mu^3 - mu * sigma^2) / sigma^2
  beta <- (mu - 2 * mu^2 + mu^3 - sigma^2 + mu * sigma^2) / sigma^2
  c(alpha = alpha, beta = beta)
}
moments_from_params <- function(alpha, beta) {
  mu <- alpha / (alpha + beta)
  sigma <- sqrt((alpha * beta) / ((alpha + beta)^2 * (alpha + beta + 1)))
  c(mu = mu, sigma = sigma)

lzachmann

library(tidyverse)
library(greta)
library(mvtnorm)
library(bayesplot)


# ---- set up ----
J <- 40  # number of groups
n_j <- sample.int(15, J, replace = TRUE)  # number of observations per group

lzachmann

lzachmann

library(tidyverse)


plus_minus <- function(x, scalar) c(x-scalar, x, x+scalar)

scenarios <- expand.grid(
  b0 = 2.86 %>% plus_minus(1),  # intercept
  b1 = -.28 %>% plus_minus(.5)  # the effect of translocation
) %>% mutate(scenario = 1:n())

lzachmann

biv_dist_plot <- function(df, x=NULL, y=NULL, bins=50) {
  library(ggplot2)
  library(gridExtra)
  

  h_top <- ggplot(data=df, aes_string(x=x)) + 
    geom_histogram(aes(y=..density..), fill = "white", color = "black", bins=bins) + 
    stat_density(colour = "blue", geom="line", size = 1, position="identity", 
                 show.legend=FALSE) +
    theme(axis.title.x = element_blank())

lzachmann

data {
  int<lower = 1> n;
  int<lower = 1> p;
  matrix[n, p] X;
  int<lower = 0> y[n];
  int<lower = 1> m;
  int<lower = 1, upper = n> n_each[m];
  int<lower = 1, upper = n> idx[m, 2]; // start and stop indices for plots
  matrix<lower = 0>[n, n] D;
}