statwonk

# ============================================================================
# Hierarchical BNN: Combining Structure + Flexibility
# ============================================================================
# Four models compared across three DGPs:
#   BHM  — group-varying intercept + slope (linear)
#   BNN  — one hidden layer, no group structure
#   HBNN — group-varying intercepts + one hidden layer (the hybrid)
#   ENS  — 50/50 mixture of BHM + BNN posterior predictive draws
#
# Three DGPs:

statwonk

# ============================================================================
# BHM vs BNN vs Ensemble: Out-of-Sample CRPS Comparison
# ============================================================================
# Three cases demonstrating the "no free lunch" principle:
#   Case 1: BHM wins when group structure is describable
#   Case 2: BNN wins when structure is complex/unstructured
#   Case 3: Their ensemble minimizes out-of-sample error across both DGPs
#
# Requirements: cmdstanr, posterior, scoringRules, dplyr, tibble, ggplot2
# ============================================================================

statwonk

# Load necessary libraries
library(gamlss)
library(rstan)
library(ggplot2)

# Step 1: Generate Johnson's SU distributed data
set.seed(42)
N <- 10000  # Sample size
gamma <- 0.5  # Shape parameter (nu)
delta <- 1    # Shape parameter (tau)

statwonk

library(tidyverse)
library(sandwich)
library(lmtest)
library(lme4)

patients <- 3e3

seq_len(patients) %>%
    map_df(function(id) {
        rnorm(1, 0, 1) -> patient_effect

statwonk

import pandas as pd
import statsmodels.api as sm
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression

class AIDataEngineer:
    def __init__(self):
        self.data = None
        self.model = None
        self.X_train = None

statwonk

library(tidyverse)
library(brms)
library(tidybayes)
# https://docs.google.com/spreadsheets/d/1vG0WdjZaYlS4_7_if-OaE3uMv3aX6zfvQMx5JvDREF0/edit?usp=sharing

# Chi.sq
prop.test(c(103, 135), c(5947, 5609),
          alternative = "less",
          conf.level = 0.99)

statwonk

library(tidyverse)
read_csv("~/Downloads/Religion by Firearm Deaths Data - Firearm Deaths by State.csv",
         skip = 4) %>%
  janitor::clean_names() -> firearm_deaths

read_csv("~/Downloads/Religion by Firearm Deaths Data - Religiousity by State.csv",
         skip = 3) %>%
  janitor::clean_names() %>%
  rename(state = region) %>%
  mutate(religious = 1 - irreligion_percent/100) %>%

statwonk

library(tidyverse)

# CPI
# https://data.bls.gov/timeseries/CUSR0000SA0&output_view=pct_1mth
readxl::read_xlsx("~/Downloads/SeriesReport-20220817182003_bb0f80.xlsx", skip = 10) %>%
  mutate(date = seq.POSIXt(as.POSIXct("1957-02-01"), by = "month", length.out = n())) %>%
  select(date, cpi = Value) %>%
  inner_join(
    # Short-term Interest Rates - https://fred.stlouisfed.org/series/DGS1MO
    read_csv("~/Downloads/DGS1MO (1).csv", col_types = cols(

statwonk

library(tidyverse)

# CPI
# https://data.bls.gov/timeseries/CUSR0000SA0&output_view=pct_1mth
readxl::read_xlsx("~/Downloads/SeriesReport-20220807120638_733b06.xlsx", skip = 10) %>%
  mutate(date = seq.POSIXt(as.POSIXct("1947-02-01"), by = "month", length.out = n())) %>%
  select(date, cpi = Value) -> cpi

cpi %>%
  filter(date >= as.POSIXct("1949-01-01")) %>%

statwonk

library(tidyverse)
library(forecast)
library(urca)
library(dynlm)
library(lmtest)
library(sandwich)

# https://data.bls.gov/timeseries/CUSR0000SA0&output_view=pct_1mth
read_csv("~/Downloads/inflation.csv", skip = 11) %>%
  janitor::clean_names() %>%