kagaya

kagaya

functions {
  real f(real mu, real s, real x){
      return(normal_lpdf(x | mu, s));
      }
     
  real log_lik_Simpson(real mu, real s, real a, real b, int M) {
    vector[M+1] lp;
    real h;
    h = (b-a)/M;
    lp[1] = f(mu, s, a);

kagaya

functions {
  real f(real mu, real s, real x){
      // separately performing numerical integration to make this computing fast
      return(normal_lpdf(x | mu, s));
      }
     
  real log_lik_Simpson(real mu, real s, real a, real b, int M) {
    // numerical integration using the Simpson's rule
    vector[M+1] lp;
    real h;

kagaya

library(rstan)
library(tidyverse)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

## preprocessing                     
d <- read_csv(file='Dromiidae.csv') # https://gist.github.com/kagaya/0d309397300b7e8294fe53c44b6fc525
d <- d %>%
  select(c('shell_width', 'choice',

kagaya

library(tidyverse)
library(ggrepel)
library(rstan)
expose_stan_functions("my_stan_functions.stan") # https://gist.github.com/kagaya/c726f16d82b80026bb1924b408a72b5c
source("utility.R") # https://gist.github.com/kagaya/60c4190ac306840daee54115b3c315c3

plot_choice_random_intercept <- function(d, fit){
  # mu[n,2] =  bias_l[ID[n]] + cwl*C_width[n] + ll*Leg_lack[n];
  # mu[n,3] =  bias_no0 + cwno*C_width[n] + lno*Leg_lack[n]; 
  ms <- rstan::extract(fit)

kagaya

waic <- function(log_likelihood) {
  # from https://gist.github.com/MatsuuraKentaro/3f6ae5863e700f5039c19e36a9bdf646
  training_error <- - mean(log(colMeans(exp(log_likelihood))))
  functional_variance_div_N <- mean(colMeans(log_likelihood^2) - colMeans(log_likelihood)^2)
  waic <- training_error + functional_variance_div_N
  return(waic)
}

kagaya

library(tidyverse)
library(grid)
library(gtable)


plot_segment_across_facets <- function(p, from=1, to=2, 
                                       from_point_id=1, 
                                       to_point_id=1, 
                                       plotout = F,
                                       gp=gpar(lty=1, alpha=0.5)){

kagaya

---
title: "an example of model comparison by WAIC (draft)"
author: "Katsushi Kagaya"
date: "`r Sys.time()`"
output:
  html_document: default
---

# motivation and aims

kagaya

functions {
  vector my_student_t_rng(int N, real nu, real mu, real sigma){
    vector[N] Y;
    for (n in 1:N){
      Y[n] = student_t_rng(nu, mu, sigma);
    }
    return(Y);
  }

  vector my_ZIP_rng(int N, real q, real lambda){

kagaya