bi_theor <- function(mu,sigma,n_stkp,betr_niveau,bi_aantal,aantal_noteren_in_tabel,toon) {
if(sigma > 0) {   
  if(betr_niveau >= 0 & betr_niveau < 1) {
    if(all(n_stkp == as.integer(n_stkp)) & n_stkp >=2) { 
      if(all(bi_aantal == as.integer(bi_aantal)) & bi_aantal >=1) { 
        if(all(aantal_noteren_in_tabel == as.integer(aantal_noteren_in_tabel)) & aantal_noteren_in_tabel >=1) { 
          if(toon %in% 1:3) {      
            library(ggplot2)
            library(dplyr)
            library(glue)
            library(gt)
            library(gtExtras)
            mu1 <- round(mu, 2)
            sigma1 <- round(sigma, 2) 
            st_err <- sigma/sqrt(n_stkp)
            st_err1 <- round(st_err, 2)
            nivperc <- trunc(betr_niveau*10^4)/10^2
            ymax <- 1/(st_err*(sqrt(2*pi)))
            ymax1 <- 1.1*ymax
            tky <- round(c(0, ymax*exp(-0.5), ymax*exp(-2), ymax) , 3)
            sim_matrix <- matrix(rnorm(n_stkp*bi_aantal, mu , sigma), ncol = bi_aantal)
            sim_mean <- apply(sim_matrix, 2, mean)
            e_margin <- qnorm((1+nivperc/100)/2)*st_err
            sim_ll <- sim_mean - e_margin
            sim_ul <- sim_mean + e_margin
            sim_nr <- 1:bi_aantal
            sim_miss <- sim_ul < mu | sim_ll > mu
            tel <- bi_aantal - sum(sim_miss)  
            sim_miss <- as.factor(sim_miss)
            perc <- round(tel/bi_aantal, 2)*100
            step <- 0.9*ymax1/(bi_aantal + 1)
            sim_id <- seq(-step, -0.9*ymax1 + step,by= -step)
            xafstbi <- max(mu - min(sim_ll) , max(sim_ul - mu)) 
            xlow <- min((mu - 3.2*st_err), mu - xafstbi)
            xhigh <- max((mu + 3.2*st_err), mu + xafstbi)  
            xlow1 <- xlow - (xhigh-xlow)*0.015
            xhigh1 <- xhigh + (xhigh-xlow)*0.015
            tkx <- round(seq(from=(mu - 6*st_err), to=(mu + 6*st_err), by=st_err) , 2)
            myfun <- function(x) {dnorm(x, mean = mu, sd = st_err)}
            df <- data.frame(sim_mean, sim_id, sim_ll, sim_ul, sim_miss, sim_nr)
            p <- ggplot(df, aes(x = sim_mean, y = sim_id, color = sim_miss)) +
              coord_cartesian(xlim =c(xlow1,xhigh1), ylim =c(-ymax1,ymax1), expand = FALSE) +
              scale_x_continuous(breaks = tkx) +
              scale_y_continuous(breaks = tky) + 
              scale_color_manual(values = c( "FALSE" = "grey65", "TRUE" = "red")) +    
              labs(title = bquote( ~ .(nivperc) ~ "% betrouwbaarheidsinterval") , subtitle = "voor het populatiegemiddelde 𝝻", y= "y", x = "x") +
              stat_function(fun = myfun , color = "black") +
              geom_point(size = 0.9) +
              geom_errorbarh(aes(xmin = sim_ll, xmax = sim_ul, y = sim_id, color = sim_miss),linewidth=0.25, height=0) +
              theme_bw() +    
              theme(
                legend.position="none" ,
                panel.grid = element_line(color = "gray95") ,
                aspect.ratio = 1 ,
                plot.title = element_text(hjust = 0.5),
                plot.subtitle = element_text(size = 13 , hjust = 0.5),
                axis.title.x = element_text(hjust = 1),
                axis.title.y = element_text(vjust = 1, angle =0)) +
              annotate("segment", x=xlow1, xend=xhigh1, y=0, yend=0, colour="black") +
              annotate("segment", x=mu, xend=mu, y=-0.9*ymax1, yend=0.025*ymax1, linewidth=0.5, lty=2, color = "grey50") +
              annotate("text", x = (5*mu + 7*xlow1)/12 , y =0.9*ymax , size = 3.5, label = "kansmodel \nsteekproefgemiddelde") +
              annotate("text", x = (mu + 2*xlow1)/3 , y =0.66*ymax , size = 3.5, label = bquote( ~ bar(X) ~ "  \u007e  N (" ~ mu ~ "," ~ frac(sigma, sqrt(n)) ~ " )")) +
              annotate("text", x = (mu + 2*xlow1)/3 , y =0.5*ymax , size = 3.5, label = bquote( ~ sigma ==.(sigma1)~" "~n==.(n_stkp))) +
              annotate("text", x = (mu + 2*xlow1)/3 , y =0.36*ymax , size = 3.5, label = bquote( "se("~ bar(X)~")" ==.(st_err1))) +   
              annotate("text", x = mu , y =0.1*ymax1 , size = 4.5, label = "𝛍" ) 
            if(bi_aantal >4){p <- p + annotate("text", x = mu , y =-0.94*ymax1 , size = 3.2, label = paste0(as.character(tel) , " van de " , as.character(bi_aantal) , " betrouwbaarheidsintervallen (", as.character(perc),  "%)  bevatten μ"))}
            fmarge <- paste("foutenmarge =", as.character(round(e_margin, 2)))
            bi_tbl <- select(df, sim_ll, sim_ul, sim_mean, sim_nr, sim_miss) |>
              head(aantal_noteren_in_tabel) |> 
              gt() |> cols_hide(sim_miss) |> 
              tab_header(title = html(glue("{nivperc} %  betrouwbaarheidsinterval voor 𝛍<br><font size=3.5>bij een normale populatie met 𝛍={mu1} en 𝞼={sigma1}</font>"))) |>
              tab_spanner(label = fmarge, columns = c(sim_ll, sim_ul) ) |> 
              tab_spanner(label = paste("n = ",as.character(n_stkp)), columns = sim_mean) |> 
              cols_label (
                sim_ll = html("ondergrens"),
                sim_ul = html("bovengrens"),     
                sim_mean = html("x&#772;"),
                sim_nr = html("volgnr.") )  |> 
              fmt_number(
                columns = c(sim_ll, sim_mean, sim_ul),
                use_seps = FALSE ,
                decimals = 2) |> 
              tab_style(style = cell_text(align = "center"),
                        locations = list( cells_column_labels(), cells_body())) |> 
              tab_options(data_row.padding = px(0.5),
                          heading.subtitle.font.size = 16) |>
              gt_highlight_rows(rows = (sim_miss == "TRUE"), fill = "red",
                                alpha = 0.3, font_weight = "normal") |> 
              gt_highlight_cols(sim_mean, fill = "lightblue",
                                alpha = 0.4, font_weight = "normal") |> 
              gt_highlight_cols(sim_nr, fill = "grey",
                                alpha = 0.05, font_weight = "normal") 
            if(bi_aantal == 1){bi_tbl <- bi_tbl |> cols_hide(sim_nr) }
            if(toon == 3){print(p)
              bi_tbl
            }else if(toon ==2) {bi_tbl 
            }else if(toon ==1) {print(p)}
          }else{cat("Voor toon <-... mag je alleen maar \nhet cijfer 1 of 2 of 3 invullen. Pas aan.")  }
        }else{cat("Tik een geheel getal groter dan 0 voor het aantal intervallen \nte noteren in de tabel: aantal_noteren_in_tabel. Pas aan.") }
      }else{cat("Tik een geheel getal groter dan 0 voor het aantal \nop te stellen betrouwbaarheidsintervallen: bi_aantal. Pas aan.") }
    }else{cat("Tik een geheel getal groter dan 1 voor \nde steekproefgrootte: n_stkp. Pas aan.") }
  }else{cat("Een betrouwbaarheidsniveau ligt tussen 0 en 1. \nBovendien is een 100 % betrouwbaarheidsniveau niet zinvol \nwant dan krijg je (-∞ , +∞) als betrouwbaarheidsinterval. \nPas dus aan zodat  0 ≤ betr_niveau < 1.") }
}else{cat("sigma moet positief zijn. Pas aan.") }
} 