

setwd("worktH_")


options(width=1000,max.print=100000,cmdstanr_max_rows=100000)

d <- read.csv("data_example.csv",fileEncoding="shift-jis")

dat <- list(
Ndata=nrow(d),
Nsp=max(d$Species_index),
Nsite=max(d$Site_index),
Sd_palatability_error=0.1,
Sd_pressure_error=0.1,
Species_index=d$Species_index,
Palatability_value_input=d$Palatability_value_input,
Palatability_value_occurence_input=d$Palatability_value_occurence_input,
Site_index=d$Site_index,
Pressure_value_input=d$Pressure_value_input,
Pressure_value_occurence_input=d$Pressure_value_occurence_input,
Examined=d$S,
Browsed=d$H
)

library(cmdstanr)

model <- cmdstan_model("palatability_cmd.stan")

fit <- model$sample(dat,iter_warmup=20000,iter_sampling=2000,thin=1, chains=16,parallel_chains=16, refresh=200)


map <- function(z){density(z)$x[which.max(density(z)$y)]}
quantile95 <- function(x){quantile(x,probs=c(0.005,0.025,0.05,0.25,0.5,0.75,0.95,0.975,0.9995),names=TRUE)}
ess_tail95 <- function(x){min(posterior::ess_quantile(x,probs=c(0.025,0.975),names=TRUE))}
r_hat <- posterior::rhat
pd_out <- c("mean","sd","map","quantile95","ess_bulk", "ess_tail95","r_hat")


sink("data_example.txt")
cat("data_example\n")
#fit$print(,pd_out)

cat("Palatability by Byesian logistic regression\n")
fit$print("palatability",pd_out)
fit$print("pressure",pd_out)
fit$print("regression_intercept",pd_out)
cat(" \n")
cat("Palatability by original scale\n")
fit$print("original_palatability",pd_out)
fit$print("original_threshold_palatability",pd_out)
cat(" \n")
cat(" \n")
cat(" \n")
print(fit)
sink()



