#########################################################################
##       R PROGRAM: Plot_LCM_priors.R
##
##         PROJECT: Evaluation of multiple-encounter population size
##                  estimators
##
##     DESCRIPTION: Plot the priors for alpha and effect of LCM K
##                  parameter
##
##      WRITTEN BY: Steve Gutreuter, CDC/CGH/DGHT Statistics, Estimation
##                                   and Modeling Team
##                  sgutreuter@cdc.gov
##
##            DATE: 2021-10-16
##
##      DISCLAIMER: Although this program has been used by the Centers
##                  for Disease Control & Prevention (CDC), no warranty,
##                  expressed or implied, is made by the CDC or the U.S.
##                  Government as to the accuracy and functioning of the
##                  program and related program material nor shall the
##                  fact of distribution constitute any such warranty,
##                  and no responsibility is assumed by the CDC in
##                  connection therewith.
##
#########################################################################
library(tidyverse)
library(ggplot2)

basepath <- file.path(Sys.getenv("PROJ"), "PSE/PSEsim")
workpath <- file.path(basepath, "R")
datapath <- file.path(basepath, "data")
outpath <- file.path(basepath, "output")
setwd(workpath)
## source(file.path(workpath, "PSE_sim_functions.R"))
## dataname <- file.path(datapath, "PSEsimSamples_alt.Rdata")
## attach(dataname)
theme_classic( )
xaxmax <- 150
alpha <- (1:xaxmax) / 10
fx2 <- dgamma(alpha, 0.025, 0.025)
fx1 <- dgamma(alpha, 0.25, 0.25)
fx3 <- dgamma(alpha, 1, 1)
labs <- c(rep("0.25, 0.25", xaxmax), rep("0.025, 0.025", xaxmax),
          rep("1.00, 1.00", xaxmax))
df1 <- data.frame(x = rep(alpha, 3), y = c(fx1, fx2, fx3), parms = labs)
ggplot(df1, aes(x, y, linetype = parms)) +
    theme_classic( ) +
    scale_linetype_manual(values = c(2,1,4)) +
    theme(axis.text.y =  element_text(size = 12),
          axis.text.x =  element_text(size = 12),
          axis.title =  element_text(size = 15),
          legend.position = c(0.5, 0.7),
          legend.title = element_text(size =  12),
          legend.text = element_text(size =  12)) +
    geom_line(aes(linetype = parms) ) +
    labs(x = expression(alpha),
         y = expression(paste("f(", alpha,  ")")),
         linetype = "Parameters:\nscale, rate" ) +
    ylim(0, 1)
ggsave(file.path(outpath, "graphs/Default_LCM_alpha_prior.png"),
       device = "png")

################################   END of FILE   ########################
