diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000..f1949a5 Binary files /dev/null and b/.DS_Store differ diff --git a/.RData b/.RData new file mode 100644 index 0000000..72fd143 Binary files /dev/null and b/.RData differ diff --git a/.gitignore b/.gitignore index 8c1d3af..b390bea 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ .Rhistory test.R .lintr -*.gz \ No newline at end of file +*.gz +.Rproj.user diff --git a/ProbDistCalc_RShiny.Rproj b/ProbDistCalc_RShiny.Rproj new file mode 100644 index 0000000..8e3c2eb --- /dev/null +++ b/ProbDistCalc_RShiny.Rproj @@ -0,0 +1,13 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX diff --git a/distribution_info.yaml b/distribution_info.yaml index 05bbe65..9f09314 100644 --- a/distribution_info.yaml +++ b/distribution_info.yaml @@ -216,9 +216,15 @@ Chi Square Non Central Distribution: Circle Distribution: id: 15 name: Circle - inputNames: CircleRadius - labels: radius - defaultValues: 2.0 + inputNames: + - CircleRadius + - CircleCenter + labels: + - radius + - center + defaultValues: + - 2.0 + - 0.0 hasImplementation: yes isWithSD: no fitFunc: fitCircle diff --git a/fitFunctions.R b/fitFunctions.R index 07c4f22..64816d2 100644 --- a/fitFunctions.R +++ b/fitFunctions.R @@ -1,27 +1,52 @@ -# FIXME: This is not working -# fitAndersonDarling <- function(dataset) { -# fitDistModel <- fitdist(dataset, "andersonDarling") -# return(fitDistModel) -# } -# FIXME: This is not working -# fitArcSine <- function(dataset) { -# fitDistModel <- fitdist(dataset, "arcsine") -# return(fitDistModel) -# } - -# FIXME: This is not working -# fitBenford <- function(dataset) { -# rounded_data <- round(dataset) -# fitDistModel <- fitdist(rounded_data, "benford") -# return(fitDistModel) -# } + fitAndersonDarling <- function(dataset) { + fitDistModel <- fitdist(dataset, "andersonDarling") + return(fitDistModel) + } + + +fitArcSine <- function(dataset) { + dataset_scaled <- (dataset - min(dataset)) / (max(dataset) - min(dataset)) + dcustom <- function(x, a, b) { + ifelse(x > a & x < b, 1 / (pi * sqrt((x - a) * (b - x))), 0) + } + + pcustom <- function(q, a, b) { + if (q > a & q < b) { + return(asin(sqrt((q - a) / (b - a))) / (pi / 2)) + } else if (q <= a) { + return(0) + } else { + return(1) + } + } + + qcustom <- function(p, a, b) { + if (p >= 0 & p <= 1) { + return(a + (b - a) * sin((pi / 2) * p)^2) + } else if (p < 0) { + return(a) + } else { + return(b) + } + } + fitDistModel <- fitdist( dataset_scaled, "custom", method = "mle", start = list(a = 0, b = 1)) + return(fitDistModel) + } + +fitBenford <- function(dataset) { + single_cell_table <- data.frame("estimate" = 1) + row.names(single_cell_table) <- "Benfn" + + return(single_cell_table) +} + +fitBernoulli <- function(dataset) { + + p_hat <- mean(dataset) + return(list(estimate = p_hat)) +} -# FIXME: This is not working -# fitBernoulli <- function(dataset) { -# fitDistModel <- fitdist(dataset, "bernoulli") -# return(fitDistModel) -# } fitBeta <- function(dataset) { @@ -43,6 +68,42 @@ fitBeta <- function(dataset) { # return(fitDistModel) # } + +library(extraDistr) +library(bbmle) + # Round to nearest integer + +fitBetaBinomial <- function(dataset) { + dataset <- round(dataset) + size <- ceiling(max(dataset)) + + # Define the negative log-likelihood function for Beta-Binomial + neg_log_likelihood <- function(alpha, beta) { + if (alpha <= 0 || beta <= 0) return(Inf) + -sum(dbbinom(dataset, size = size, alpha = alpha, beta = beta, log = TRUE)) + } + + # Fit the model using MLE + fit <- mle2( + neg_log_likelihood, + start = list(alpha = 1, beta = 1), + method = "L-BFGS-B", + lower = c(0.0001, 0.0001), + control = list(maxit = 1000) + ) + fit_df <- data.frame( + estimate = c( size,coef(fit)["alpha"], coef(fit)["beta"]), + row.names = c( "n", "alpha", "beta") + ) + + return(fit_df) +} + + + + + + # FIXME: This is not working # fitBinomial <- function(dataset) { # rounded_data <- round(dataset) @@ -85,11 +146,38 @@ fitCauchy <- function(dataset) { # return(fitDistModel) # } -# FIXME: This is not working -# fitCircle <- function(dataset) { -# fitDistModel <- fitdist(dataset, "circle") -# return(fitDistModel) -# } +library(minpack.lm) + +fitCircle <- function(dataset) { + # Define residual function to minimize + circleResiduals <- function(params) { + r <- params[1] # radius of the circle + c <- params[2] # center of the circle + abs(dataset - c) - r # residuals: difference from each point to radius centered at c + } + + # Provide initial guesses for radius and center + start_center <- mean(dataset) # initial guess for center + start_radius <- max(abs(dataset - start_center)) # initial guess for radius + + # Fit radius and center using nonlinear least squares + fitResult <- nls.lm(par = c(start_radius, start_center), fn = circleResiduals, control = list(maxiter = 100)) + + # Extract the fitted radius and center + fitted_radius <- fitResult$par[1] + fitted_center <- fitResult$par[2] + + + # Return a data frame with one row and one column containing the fitted radius + fit_df <- data.frame( + estimate = c( fitted_radius, fitted_center), + row.names = c( "radius", "center") + ) + return(fit_df) + + + + } # Not tested fitContinuousUniform <- function(dataset) { @@ -97,28 +185,194 @@ fitContinuousUniform <- function(dataset) { return(fitDistModel) } -# Not tested +coupon_likelihood <- function(params, data) { + N <- params[1] # Total population size (number of distinct items) + k <- params[2] # Number of distinct items needed + + # Calculate expected trials based on harmonic sum approximation + expected_trials <- N * sum(1 / (N - (0:(k - 1)))) + + # Round data to integer values for Poisson + rounded_data <- round(data) + + # Calculate the negative log-likelihood using the rounded data + -sum(dpois(rounded_data, lambda = expected_trials, log = TRUE)) +} +# Not tested +coupon_likelihood <- function(params, data) { + N <- params[1] # Total population size (number of distinct items) + k <- params[2] # Number of distinct items needed + + # Calculate expected trials based on harmonic sum approximation + expected_trials <- N * sum(1 / (N - (0:(k - 1)))) + + # Round data to integer values for Poisson + rounded_data <- round(data) + + # Calculate the negative log-likelihood using the rounded data + -sum(dpois(rounded_data, lambda = expected_trials, log = TRUE)) +} fitCoupon <- function(dataset) { - fitDistModel <- fitdist(dataset, "coupon") - return(fitDistModel) + # Initial guesses for N and k + start_params <- c(N = max(dataset), k = round(mean(dataset) / 2)) + + # Perform optimization with integer constraint on k + fit_result <- optim( + par = start_params, + fn = coupon_likelihood, + data = dataset, + method = "L-BFGS-B", + lower = c(1, 1), + upper = c(Inf, Inf) + ) + + # Extract fitted parameters, rounding k to the nearest integer + N_estimate <- round(fit_result$par[1]) + k_estimate <- round(fit_result$par[2]) # Ensure integer k + + # Return the results as a data frame + fit_df <- data.frame( + estimate = c(N_estimate, k_estimate), + row.names = c("Population Size", "Number of distinct values needed") + ) + + return(fit_df) } # Not tested fitDie <- function(dataset) { - fitDistModel <- fitdist(dataset, "die") - return(fitDistModel) + getDensity_Die <- function(x) { + k <- round(x) + if (n_Die == 0) { + # FAIR Die + if (k < 1 || k > 6) { + return(0) + } else { + return(1/6) + } + } else if (n_Die == 1) { + # FLAT16 Die + if (k < 1 || k > 6) { + return(0) + } else if (k == 1 || k == 6) { + return(1/4) + } else { + return(1/8) + } + } else if (n_Die == 2) { + # FLAT25 Die + if (k < 1 || k > 6) { + return(0) + } else if (k == 2 || k == 5) { + return(1/4) + } else { + return(1/8) + } + } else if (n_Die == 3) { + # FLAT34 Die + if (k < 1 || k > 6) { + return(0) + } else if (k == 3 || k == 4) { + return(1/4) + } else { + return(1/8) + } + } else if (n_Die == 4) { + # LEFT Die + if (k < 1 || k > 6) { + return(0) + } else if (k == 1) { + return(1/21) + } else if (k == 2) { + return(2/21) + } else if (k == 3) { + return(3/21) + } else if (k == 4) { + return(4/21) + } else if (k == 5) { + return(5/21) + } else { + return(6/21) + } + } else if (n_Die == 5) { + # RIGHT Die + if (k < 1 || k > 6) { + return(0) + } else if (k == 1) { + return(6/21) + } else if (k == 2) { + return(5/21) + } else if (k == 3) { + return(4/21) + } else if (k == 4) { + return(3/21) + } else if (k == 5) { + return(2/21) + } else { + return(1/21) + } + } + } + best_die_type <- NULL + lowest_residual <- Inf + residuals_list <- list() + + + for (die_type in 0:5) { + n_Die <<- die_type + expected_density <- sapply(dataset, getDensity_Die) + residuals <- (1 / length(dataset) - expected_density)^2 + residual_sum <- sum(residuals) + residuals_list[[paste0("Die_Type_", die_type)]] <- residual_sum + + if (residual_sum < lowest_residual) { + lowest_residual <- residual_sum + best_die_type <- die_type + } + } + fit_df <- data.frame( + estimate = best_die_type, + row.names = 'n (0 <= n <= 5)' + ) + return(fit_df) } # Not tested fitDiscreteArcSine <- function(dataset) { - fitDistModel <- fitdist(dataset, "disarc") + dataset <- round(dataset) + darcsine <- function(x, a, b) { + if (a >= b) stop("Invalid range: a must be less than b") + ifelse(x >= a & x <= b, 1 / (pi * sqrt((x - a) * (b - x))), 0) + } + + # CDF of the arcsine distribution + parcsine <- function(q, a, b) { + if (a >= b) stop("Invalid range: a must be less than b") + ifelse( + q < a, 0, + ifelse( + q > b, 1, + (2 / pi) * asin(sqrt((q - a) / (b - a))) + ) + ) + } + + # Quantile function (inverse CDF) of the arcsine distribution + qarcsine <- function(p, a, b) { + if (a >= b) stop("Invalid range: a must be less than b") + if (p < 0 | p > 1) stop("Probability p must be in [0, 1]") + a + (b - a) * sin((p * pi) / 2)^2 + } + + fitDistModel <- fitdist( dataset, "arcsine", method = "mle", start = list(a = 1, b = 2)) return(fitDistModel) } # Not tested fitDiscreteUniform <- function(dataset) { - rounded_data <- round(dataset) - fitDistModel <- fitdist(rounded_data, "discunif") + #rounded_data <- round(dataset) + fitDistModel <-fitdist(dataset, "unif", method = c("mme"), + start=NULL, fix.arg=NULL, discrete = TRUE) return(fitDistModel) } @@ -128,6 +382,23 @@ fitDiscreteUniform <- function(dataset) { # return(fitDistModel) # } +fitErlang <- function(dataset) { + mean_data <- mean(dataset) + var_data <- var(dataset) + + # Rough estimate of the shape parameter + myshape <- round(mean_data^2/var_data) + + + fit_df <- data.frame( + estimate = c( mean_data/myshape, myshape), + row.names = c( "scale", "shape") + ) + return(fit_df) + +} + + # FIXME: This is not working # fitError <- function(dataset) { # fitDistModel <- fitdist(dataset, "error") diff --git a/plotlyFunctions/Circle.R b/plotlyFunctions/Circle.R index 9df32bf..82dea4d 100644 --- a/plotlyFunctions/Circle.R +++ b/plotlyFunctions/Circle.R @@ -1,68 +1,69 @@ -dCircle <- function(x, radius) { - ifelse((-radius <= x & x <= radius), (2 * sqrt(radius * radius - x * x) / (pi * radius * radius)), 0) +dCircle <- function(x, radius, center) { + ifelse((-radius + center <= x & x <= radius + center), + (2 * sqrt(radius * radius - (x - center) * (x - center)) / (pi * radius * radius)), + 0) } -pCircle <- function(x, radius) { - 0.5 + asin(x / radius) / pi + x * sqrt(1 - x * x / (radius * radius)) / (pi * radius) +pCircle <- function(x, radius, center) { + adjusted_x <- x - center + 0.5 + asin(adjusted_x / radius) / pi + adjusted_x * sqrt(1 - adjusted_x * adjusted_x / (radius * radius)) / (pi * radius) } -plotlyCircleDistribution <- function(plotrange, input, distType, probrange) { - xseq <- seq( - min(0, as.numeric(plotrange[1])), max(as.numeric(plotrange[2]), 10), - 0.01 - ) - f15 <- 0 - if (input$FunctionType == "PDF/PMF") { - f15 <- dCircle(xseq, as.numeric(input$CircleRadius)) - graphtype <- "PMF" - } else if (input$FunctionType == "CDF/CMF") { - f15 <- pCircle(xseq, as.numeric(input$CircleRadius)) - graphtype <- "CMF" - } else { - graphtype <- "" +plotlyCircleDistribution <- function(plotrange, input, distType, probrange, center) { + xseq <- seq( + min(0, as.numeric(plotrange[1])), max(as.numeric(plotrange[2]), 10), + 0.01 + ) + f15 <- 0 + if (input$FunctionType == "PDF/PMF") { + f15 <- dCircle(xseq, as.numeric(input$CircleRadius), as.numeric(input$CircleCenter)) + graphtype <- "PMF" + } else if (input$FunctionType == "CDF/CMF") { + f15 <- pCircle(xseq, as.numeric(input$CircleRadius), as.numeric(input$CircleCenter)) + graphtype <- "CMF" + } else { + graphtype <- "" + } + if (graphtype != "") { + xsize <- length(xseq) + colors <- c(rep("rgb(31, 119, 180)", xsize)) + for (index in 1:xsize) { + if (xseq[index] >= round(probrange[1], 0) && xseq[index] <= round(probrange[2], 0)) { + colors[index] <- "rgb(255, 127, 14)" + } } - if (graphtype != "") { - xsize <- length(xseq) - colors <- c(rep("rgb(31, 119, 180)", xsize)) - for (index in 1:xsize) { - if (xseq[index] >= round(probrange[1], 0) && xseq[index] <= round( - probrange[2], - 0 - )) { - colors[index] <- "rgb(255, 127, 14)" - } - } - - prob <- pCircle(as.numeric(probrange[2]), as.numeric(input$CircleRadius)) - pCircle(as.numeric(probrange[1]), as.numeric(input$CircleRadius)) - - fig <- plot_ly( - x = xseq, y = f15, name = distType, type = "scatter", mode = "lines", - hoverinfo = "xy" - ) - - fig <- fig %>% - add_trace( - x = xseq, y = f15, name = paste("Probability = ", prob, sep = ""), - hoverinfo = "name", fill = "tozeroy", fillcolor = "rgba(255, 212, 96, 0.5)" - ) - fig <- fig %>% - plotly::layout( - title = paste(distributions[15], " - ", graphtype, sep = ""), - hovermode = "x", hoverlabel = list(namelength = 100), yaxis = list( - fixedrange = TRUE, - zeroline = TRUE, range = c(min(f15), max(f15)), type = "linear" - ), - xaxis = list( - showticklabels = TRUE, title = "* All x values rounded to nearest integers", - zeroline = TRUE, showline = TRUE, showgrid = TRUE, linecolor = "rgb(204, 204, 204)", - linewidth = 2, mirror = TRUE, fixedrange = TRUE, range = c( - plotrange[1], - plotrange[2] - ) - ), showlegend = FALSE - ) - fig <- fig %>% - config(editable = FALSE) - fig - } -} + + prob <- pCircle(as.numeric(probrange[2]), as.numeric(input$CircleRadius), as.numeric(input$CircleCenter)) - + pCircle(as.numeric(probrange[1]), as.numeric(input$CircleRadius), as.numeric(input$CircleCenter)) + + fig <- plot_ly( + x = xseq, y = f15, name = distType, type = "scatter", mode = "lines", + hoverinfo = "xy" + ) + + fig <- fig %>% + add_trace( + x = xseq, y = f15, name = paste("Probability = ", prob, sep = ""), + hoverinfo = "name", fill = "tozeroy", fillcolor = "rgba(255, 212, 96, 0.5)" + ) + fig <- fig %>% + plotly::layout( + title = paste(distributions[15], " - ", graphtype, sep = ""), + hovermode = "x", hoverlabel = list(namelength = 100), yaxis = list( + fixedrange = TRUE, + zeroline = TRUE, range = c(min(f15), max(f15)), type = "linear" + ), + xaxis = list( + showticklabels = TRUE, title = "* All x values rounded to nearest integers", + zeroline = TRUE, showline = TRUE, showgrid = TRUE, linecolor = "rgb(204, 204, 204)", + linewidth = 2, mirror = TRUE, fixedrange = TRUE, range = c( + plotrange[1], + plotrange[2] + ) + ), showlegend = FALSE + ) + fig <- fig %>% + config(editable = FALSE) + fig + } +} \ No newline at end of file diff --git a/renv.lock b/renv.lock index 7d826f2..c7ffd11 100644 --- a/renv.lock +++ b/renv.lock @@ -66,7 +66,7 @@ }, "DescTools": { "Package": "DescTools", - "Version": "0.99.50", + "Version": "0.99.58", "Source": "Repository", "Repository": "RSPM", "Requirements": [ @@ -167,7 +167,7 @@ }, "Matrix": { "Package": "Matrix", - "Version": "1.6-0", + "Version": "1.7-0", "Source": "Repository", "Repository": "CRAN", "Requirements": [ @@ -215,7 +215,7 @@ }, "RcppArmadillo": { "Package": "RcppArmadillo", - "Version": "0.12.6.4.0", + "Version": "14.2.2-1", "Source": "Repository", "Repository": "RSPM", "Requirements": [ @@ -229,7 +229,7 @@ }, "RcppEigen": { "Package": "RcppEigen", - "Version": "0.3.3.9.4", + "Version": "0.3.4.0.2", "Source": "Repository", "Repository": "CRAN", "Requirements": [ @@ -268,7 +268,7 @@ }, "VGAM": { "Package": "VGAM", - "Version": "1.1-9", + "Version": "1.1-12", "Source": "Repository", "Repository": "RSPM", "Requirements": [ @@ -475,7 +475,7 @@ }, "circular": { "Package": "circular", - "Version": "0.5-0", + "Version": "0.5-1", "Source": "Repository", "Repository": "RSPM", "Requirements": [ @@ -600,7 +600,7 @@ }, "deSolve": { "Package": "deSolve", - "Version": "1.38", + "Version": "1.40", "Source": "Repository", "Repository": "RSPM", "Requirements": [ @@ -727,7 +727,7 @@ }, "expm": { "Package": "expm", - "Version": "0.999-7", + "Version": "1.0-0", "Source": "Repository", "Repository": "RSPM", "Requirements": [ @@ -1163,7 +1163,7 @@ }, "lmom": { "Package": "lmom", - "Version": "3.0", + "Version": "3.2", "Source": "Repository", "Repository": "RSPM", "Requirements": [ @@ -1209,7 +1209,7 @@ }, "mgcv": { "Package": "mgcv", - "Version": "1.9-0", + "Version": "1.9-1", "Source": "Repository", "Repository": "CRAN", "Requirements": [ @@ -1283,7 +1283,7 @@ }, "mvtnorm": { "Package": "mvtnorm", - "Version": "1.2-3", + "Version": "1.3-2", "Source": "Repository", "Repository": "RSPM", "Requirements": [ @@ -1294,7 +1294,7 @@ }, "nlme": { "Package": "nlme", - "Version": "3.1-162", + "Version": "3.1-166", "Source": "Repository", "Repository": "CRAN", "Requirements": [ @@ -1308,7 +1308,7 @@ }, "nloptr": { "Package": "nloptr", - "Version": "2.0.3", + "Version": "2.1.1", "Source": "Repository", "Repository": "RSPM", "Requirements": [ @@ -1664,7 +1664,7 @@ }, "rgl": { "Package": "rgl", - "Version": "1.2.1", + "Version": "1.3.1", "Source": "Repository", "Repository": "CRAN", "Requirements": [ @@ -1744,7 +1744,7 @@ }, "rstpm2": { "Package": "rstpm2", - "Version": "1.6.2", + "Version": "1.6.6", "Source": "Repository", "Repository": "RSPM", "Requirements": [ @@ -2045,7 +2045,7 @@ }, "tolerance": { "Package": "tolerance", - "Version": "2.0.0", + "Version": "3.0.0", "Source": "Repository", "Repository": "CRAN", "Requirements": [ @@ -2198,7 +2198,7 @@ }, "xml2": { "Package": "xml2", - "Version": "1.3.5", + "Version": "1.3.6", "Source": "Repository", "Repository": "CRAN", "Requirements": [ diff --git a/renv/.DS_Store b/renv/.DS_Store new file mode 100644 index 0000000..cb529a3 Binary files /dev/null and b/renv/.DS_Store differ diff --git a/server.R b/server.R index bd70f09..f469497 100644 --- a/server.R +++ b/server.R @@ -1,267 +1,283 @@ -# SOCR Probability Distribution Calculator -# Version 0.9 -# Updated December 8th by Bole Li and Joonseop Kim at the University of Michigan -SOCR -# Orginally created by Jared(Tianyi) Chai - -# This is a SOCR Interactive Graphical Probability Distribution Calculator -# You can run the application by clicking -# the 'Run App' button above. - -# ----------------------- Server.R ----------------------- # -# For backend calculations -library(xml2) -library(shinyjs) -library(flexsurv) -library(vcdExtra) -library(evd) -library(DescTools) -library(shiny) -library(triangle) -library(plotly) -library(stringr) -library(VGAM) -library(BayesTools) -library(extraDistr) -library(statmod) -library(truncnorm) -library(tolerance) -library(chi) -library(Rlab) -library(shinyWidgets) -library(circular) -library(mnormt) -library(ExtDist) -library(VaRES) -library(shinyjs) -shinyjs::useShinyjs() -source("renderMainPlot.R") -source("renderProbability.R") - -shinyServer( - function(input, output, session) { - datasetShown <- reactiveVal(iris) - # ----------------------- Update Distribution Type and Function Type according to URL handle ----------------------- # - observe({ - query <- parseQueryString(session$clientData$url_search) - if (!is.null(query[["d"]])) { - updateSelectInput(session, "Distribution", selected = distributions[as.numeric(query[["d"]])]) - updateSelectInput(session, "FunctionType", selected = "PDF/PMF") - if (!is.null(query[["t"]])) { - updateSelectInput(session, "FunctionType", selected = query[["t"]]) - } - } - # ----------------------- Update Range of Probability Calculation according to Range of X ----------------------- # - updateSliderInput(session, - "probrange", - value = 0, - min = input$plotrange[1], - max = input$plotrange[2], - step = 0.01 - ) - updateNumericInput(session, - "probrangeNumMin", - value = 0, - min = input$plotrangeNumMin, - max = input$plotrangeNumMax - ) - updateNumericInput(session, - "probrangeNumMax", - value = 0, - min = input$plotrangeNumMin, - max = input$plotrangeNumMax - ) - }) - - observeEvent(input$fitParams, { - updateTextInput(session, "FunctionType", value = "PDF/PMF") - distributionInfo <- distributionInfoList[[input$Distribution]] - if (is.null(dataset)) { - showNotification("Dataset is not specified.", type = "error", duration = 2) - } else if (is.null((distributionInfo$fitFunc))) { - showNotification("Fitting this distribution is not supported yet.", type = "error", duration = 2) - } else { - fit_result <- distributionInfo$fitFunc(dataset[, input$outcome]) - for (i in 1:length(fit_result$estimate)) { - inputName <- distributionInfo$inputNames[[i]] - fitted_parameter <- round(fit_result$estimate[[i]], digits = 4) - updateTextInput(session, inputName, value = fitted_parameter) - session$sendCustomMessage("highlightTextInput", inputName) - } - if (distributionInfo$isWithSD) { - # update the plot range, make it centered at the mean - oldPlotRange <- input$plotrange - halfLength <- (oldPlotRange[2] - oldPlotRange[1]) / 2 - updateSliderInput(session, - "plotrange", - label = NULL, - value = c(fit_result$estimate[[1]] - halfLength, fit_result$estimate[[1]] + halfLength), - min = -1000, - max = 1000, - step = NULL, - timeFormat = NULL, - timezone = NULL - ) - } - } - }) - - observe({ - if (input$numericalValues == 0 && input$Distribution %in% distWithSD) { - shinyjs::enable("SDNum") - # FIXME: This is not working - # shinyjs::toggle("SDNumColumn", condition = TRUE) - } else { - shinyjs::disable("SDNum") - # FIXME: This is not working - # shinyjs::toggle("SDNumColumn", condition = FALSE) - } - }) - - # Generate text to display current parameters and their values - output$currentParameters <- renderUI( - HTML({ - distributionInfo <- distributionInfoList[[input$Distribution]] - paramValues <- lapply(seq_along(distributionInfo$labels), function(i) { - label <- distributionInfo$labels[[i]] - inputName <- distributionInfo$inputNames[[i]] - paste("", label, ": ", input[[inputName]], "") - }) - # Combine the distribution name and parameter values into a single string - paste("Current Parameters:
", paste(paramValues, collapse = "
")) - }) - ) - - observeEvent(input$CalcModelerTabsetPanel, { - if (input$CalcModelerTabsetPanel == "Modeler") { - updateTextInput(session, "FunctionType", value = "") - } - }) - - # Reactive function to read uploaded file and update dataset - observeEvent(input$file, { - req(input$file) - dataset <<- read.csv(input$file$datapath) - datasetShown(dataset) - # Update choices for selectInput widgets - updateSelectInput(session, "outcome", choices = namedListOfFeatures(), selected = NULL) - updateSelectInput(session, "indepvar", choices = namedListOfFeatures(), selected = NULL) - }) - - # ----------------------- HelpMe ----------------------- # - observeEvent(input$vh.readme, { - showModal(modalDialog( - title = "Help / ReadMe", - HTML("
- SOCR Interactive Probability Distribution Calculator [Version: V.0.9] - The SOCR RShiny probability distribution calculators provide interactive vizualizations of probability densities, - mass functions, and cumulative distributions, e.g., bivariate normal distribution. -

- Acknowledgments -

- This work is supported in part by NIH grants P20 NR015331, UL1TR002240, P30 DK089503, UL1TR002240, - and NSF grants 1916425, 1734853, 1636840, 1416953, 0716055 and 1023115. Students, trainees, scholars, - and researchers from SOCR, BDDS, MIDAS, MICHR, and the broad R-statistical computing community have contributed ideas, - code, and support. -
- Developers
- Jared (Tianyi) Chai (chtianyi@umich.edu)
- Shihang Li (shihangl@umich.edu)
- Yongxiang Zhao (zyxleo@umich.edu)
- Bole Li (boleli@umich.edu)
- Joonseop Kim (joonkim@umich.edu)
- Ivo Dinov (dinov@med.umich.edu).

-
- "), - easyClose = TRUE - )) - }) - # ----------------------- Render Metadata Information from xml Database ----------------------- # - output$MetaData <- renderPrint({ - distType <- input$Distribution - distType <- tolower(str_replace_all(distType, "[^[:alnum:]]", "")) - counter <- 0 - for (i in 1:xml_len) { - j <- 1 - while (distributions_meta[[j, i * 2 - 1]] == "name") { - if (tolower(str_replace_all(distributions_meta[[1, i * 2]], "[^[:alnum:]]", "")) == distType) { - counter <- i - break - } else { - j <- j + 1 - } - } - } - outputstring <- "" - if (counter != 0) { - row <- 1 - while (distributions_meta[[row, counter * 2 - 1]] != "" && row < xml_wid) { - outputstring <- paste(outputstring, "", distributions_meta[[row, counter * 2 - 1]], ": ", distributions_meta[[row, counter * 2]], "\n", sep = "") - row <- row + 1 - } - } - withMathJax(helpText(HTML(outputstring))) - }) - # ----------------------- Render Main Plot ----------------------- # - renderMainPlot(input, output, session, dataset) - # ----------------------- Render Implementing Message ----------------------- # - output$Implementing <- renderText({ - if (input$Distribution %in% distToImpl) { - paste("The ", input$Distribution, " is still being implemented.", sep = "") - } - }) - # ----------------------- Calculate and Render Probability ----------------------- # - renderProbability(input, output, session) - - - # Imputation of categorical variables using Mode - getmode <- function(v) { - v <- v[nchar(as.character(v)) > 0] - uniqv <- unique(v) - uniqv[which.max(tabulate(match(v, uniqv)))] - } - - # Render the DataTable dynamically based on the reactive dataset - output$tbl <- DT::renderDataTable({ - # Use isolate to prevent invalidation of the reactive expression on initial render - DT::datatable(datasetShown(), options = list(lengthChange = FALSE)) - }) - - # Regression output - output$summary <- renderPrint({ - fit <- lm(unlist(dataset[, input$outcome]) ~ unlist(dataset[, input$indepvar])) - names(fit$coefficients) <- c("Intercept", input$var2) - fitCurrent <- fit - summary(fit) - }) - - # Scatterplot output - output$scatterplot <- renderPlotly({ - plot_ly( - x = ~ unlist(dataset[, input$indepvar]), y = ~ unlist(dataset[, input$outcome]), - type = "scatter", mode = "markers", name = "Data" - ) %>% - add_lines( - x = ~ unlist(dataset[, input$indepvar]), - y = ~ (lm(unlist(dataset[, input$outcome]) ~ unlist(dataset[, input$indepvar]))$fitted.values), - mode = "lines", name = "Linear Model" - ) %>% - add_lines( - x = ~ lowess(unlist(dataset[, input$indepvar]), unlist(dataset[, input$outcome]))$x, - y = ~ lowess(unlist(dataset[, input$indepvar]), unlist(dataset[, input$outcome]))$y, - mode = "lines", name = "LOESS" - ) %>% - add_markers( - x = mean(unlist(dataset[, input$indepvar])), y = mean(unlist(dataset[, input$outcome])), - name = "Center Point", marker = list(size = 20, color = "green", line = list(color = "yellow", width = 2)) - ) %>% - layout( - title = paste0( - "lm(", input$outcome, " ~ ", input$indepvar, - "), Cor(", input$indepvar, ",", input$outcome, ") = ", - round(cor(unlist(dataset[, input$indepvar]), unlist(dataset[, input$outcome])), 3) - ), - xaxis = list(title = input$indepvar), yaxis = list(title = input$outcome) - ) - }) - } -) +# SOCR Probability Distribution Calculator +# Version 0.9 +# Updated December 8th by Bole Li and Joonseop Kim at the University of Michigan -SOCR +# Orginally created by Jared(Tianyi) Chai + +# This is a SOCR Interactive Graphical Probability Distribution Calculator +# You can run the application by clicking +# the 'Run App' button above. + +# ----------------------- Server.R ----------------------- # +# For backend calculations +library(xml2) +library(shinyjs) +library(flexsurv) +library(vcdExtra) +library(evd) +library(DescTools) +library(shiny) +library(triangle) +library(plotly) +library(stringr) +library(VGAM) +library(BayesTools) +library(extraDistr) +library(statmod) +library(truncnorm) +library(tolerance) +library(chi) +library(Rlab) +library(shinyWidgets) +library(circular) +library(mnormt) +library(ExtDist) +library(VaRES) +library(shinyjs) +library(benford.analysis) +shinyjs::useShinyjs() +source("renderMainPlot.R") +source("renderProbability.R") + +shinyServer( + function(input, output, session) { + datasetShown <- reactiveVal(iris) + # ----------------------- Update Distribution Type and Function Type according to URL handle ----------------------- # + observe({ + query <- parseQueryString(session$clientData$url_search) + if (!is.null(query[["d"]])) { + updateSelectInput(session, "Distribution", selected = distributions[as.numeric(query[["d"]])]) + updateSelectInput(session, "FunctionType", selected = "PDF/PMF") + if (!is.null(query[["t"]])) { + updateSelectInput(session, "FunctionType", selected = query[["t"]]) + } + } + # ----------------------- Update Range of Probability Calculation according to Range of X ----------------------- # + updateSliderInput(session, + "probrange", + value = 0, + min = input$plotrange[1], + max = input$plotrange[2], + step = 0.01 + ) + updateNumericInput(session, + "probrangeNumMin", + value = 0, + min = input$plotrangeNumMin, + max = input$plotrangeNumMax + ) + updateNumericInput(session, + "probrangeNumMax", + value = 0, + min = input$plotrangeNumMin, + max = input$plotrangeNumMax + ) + }) + + observeEvent(input$fitParams, { + updateTextInput(session, "FunctionType", value = "PDF/PMF") + distributionInfo <- distributionInfoList[[input$Distribution]] + if (input$Distribution == "Bernoulli Distribution") { + if (is.null(dataset)) { + showNotification("Dataset is not specified.", type = "error", duration = 2) + } else if (!all(dataset[, input$outcome] %in% c(0, 1))) { + # Show a warning if the dataset is not binary + showNotification("Dataset must be binary for the Bernoulli distribution.", type = "error", duration = 4) + } + } + # Check if fitting is supported for other distributions + + + if (is.null(dataset)) { + showNotification("Dataset is not specified.", type = "error", duration = 2) + } else if (is.null((distributionInfo$fitFunc))) { + showNotification("Fitting this distribution is not supported yet.", type = "error", duration = 2) + } else { + fit_result <- distributionInfo$fitFunc(dataset[, input$outcome]) + #if (input$Distribution == "Circle Distribution") { + # selected_data <- dataset[, c(input$outcome, input$indepvar)] + # fit_result <- distributionInfo$fitFunc(selected_data) + #} + for (i in 1:length(fit_result$estimate)) { + inputName <- distributionInfo$inputNames[[i]] + fitted_parameter <- round(fit_result$estimate[[i]], digits = 4) + updateTextInput(session, inputName, value = fitted_parameter) + session$sendCustomMessage("highlightTextInput", inputName) + } + if (distributionInfo$isWithSD) { + # update the plot range, make it centered at the mean + oldPlotRange <- input$plotrange + halfLength <- (oldPlotRange[2] - oldPlotRange[1]) / 2 + updateSliderInput(session, + "plotrange", + label = NULL, + value = c(fit_result$estimate[[1]] - halfLength, fit_result$estimate[[1]] + halfLength), + min = -1000, + max = 1000, + step = NULL, + timeFormat = NULL, + timezone = NULL + ) + } + } + }) + + observe({ + if (input$numericalValues == 0 && input$Distribution %in% distWithSD) { + shinyjs::enable("SDNum") + # FIXME: This is not working + # shinyjs::toggle("SDNumColumn", condition = TRUE) + } else { + shinyjs::disable("SDNum") + # FIXME: This is not working + # shinyjs::toggle("SDNumColumn", condition = FALSE) + } + }) + + # Generate text to display current parameters and their values + output$currentParameters <- renderUI( + HTML({ + distributionInfo <- distributionInfoList[[input$Distribution]] + paramValues <- lapply(seq_along(distributionInfo$labels), function(i) { + label <- distributionInfo$labels[[i]] + inputName <- distributionInfo$inputNames[[i]] + paste("", label, ": ", input[[inputName]], "") + }) + # Combine the distribution name and parameter values into a single string + paste("Current Parameters:
", paste(paramValues, collapse = "
")) + }) + ) + + observeEvent(input$CalcModelerTabsetPanel, { + if (input$CalcModelerTabsetPanel == "Modeler") { + updateTextInput(session, "FunctionType", value = "") + } + }) + + # Reactive function to read uploaded file and update dataset + observeEvent(input$file, { + req(input$file) + dataset <<- read.csv(input$file$datapath) + datasetShown(dataset) + # Update choices for selectInput widgets + updateSelectInput(session, "outcome", choices = namedListOfFeatures(), selected = NULL) + updateSelectInput(session, "indepvar", choices = namedListOfFeatures(), selected = NULL) + }) + + # ----------------------- HelpMe ----------------------- # + observeEvent(input$vh.readme, { + showModal(modalDialog( + title = "Help / ReadMe", + HTML("
+ SOCR Interactive Probability Distribution Calculator [Version: V.0.9] + The SOCR RShiny probability distribution calculators provide interactive vizualizations of probability densities, + mass functions, and cumulative distributions, e.g., bivariate normal distribution. +

+ Acknowledgments +

+ This work is supported in part by NIH grants P20 NR015331, UL1TR002240, P30 DK089503, UL1TR002240, + and NSF grants 1916425, 1734853, 1636840, 1416953, 0716055 and 1023115. Students, trainees, scholars, + and researchers from SOCR, BDDS, MIDAS, MICHR, and the broad R-statistical computing community have contributed ideas, + code, and support. +
+ Developers
+ Jared (Tianyi) Chai (chtianyi@umich.edu)
+ Shihang Li (shihangl@umich.edu)
+ Yongxiang Zhao (zyxleo@umich.edu)
+ Bole Li (boleli@umich.edu)
+ Joonseop Kim (joonkim@umich.edu)
+ Ivo Dinov (dinov@med.umich.edu).

+
+ "), + easyClose = TRUE + )) + }) + # ----------------------- Render Metadata Information from xml Database ----------------------- # + output$MetaData <- renderPrint({ + distType <- input$Distribution + distType <- tolower(str_replace_all(distType, "[^[:alnum:]]", "")) + counter <- 0 + for (i in 1:xml_len) { + j <- 1 + while (distributions_meta[[j, i * 2 - 1]] == "name") { + if (tolower(str_replace_all(distributions_meta[[1, i * 2]], "[^[:alnum:]]", "")) == distType) { + counter <- i + break + } else { + j <- j + 1 + } + } + } + outputstring <- "" + if (counter != 0) { + row <- 1 + while (distributions_meta[[row, counter * 2 - 1]] != "" && row < xml_wid) { + outputstring <- paste(outputstring, "", distributions_meta[[row, counter * 2 - 1]], ": ", distributions_meta[[row, counter * 2]], "\n", sep = "") + row <- row + 1 + } + } + withMathJax(helpText(HTML(outputstring))) + }) + # ----------------------- Render Main Plot ----------------------- # + renderMainPlot(input, output, session, dataset) + # ----------------------- Render Implementing Message ----------------------- # + output$Implementing <- renderText({ + if (input$Distribution %in% distToImpl) { + paste("The ", input$Distribution, " is still being implemented.", sep = "") + } + }) + # ----------------------- Calculate and Render Probability ----------------------- # + renderProbability(input, output, session) + + + # Imputation of categorical variables using Mode + getmode <- function(v) { + v <- v[nchar(as.character(v)) > 0] + uniqv <- unique(v) + uniqv[which.max(tabulate(match(v, uniqv)))] + } + + # Render the DataTable dynamically based on the reactive dataset + output$tbl <- DT::renderDataTable({ + # Use isolate to prevent invalidation of the reactive expression on initial render + DT::datatable(datasetShown(), options = list(lengthChange = FALSE)) + }) + + # Regression output + output$summary <- renderPrint({ + fit <- lm(unlist(dataset[, input$outcome]) ~ unlist(dataset[, input$indepvar])) + names(fit$coefficients) <- c("Intercept", input$var2) + fitCurrent <- fit + summary(fit) + }) + + # Scatterplot output + output$scatterplot <- renderPlotly({ + plot_ly( + x = ~ unlist(dataset[, input$indepvar]), y = ~ unlist(dataset[, input$outcome]), + type = "scatter", mode = "markers", name = "Data" + ) %>% + add_lines( + x = ~ unlist(dataset[, input$indepvar]), + y = ~ (lm(unlist(dataset[, input$outcome]) ~ unlist(dataset[, input$indepvar]))$fitted.values), + mode = "lines", name = "Linear Model" + ) %>% + add_lines( + x = ~ lowess(unlist(dataset[, input$indepvar]), unlist(dataset[, input$outcome]))$x, + y = ~ lowess(unlist(dataset[, input$indepvar]), unlist(dataset[, input$outcome]))$y, + mode = "lines", name = "LOESS" + ) %>% + add_markers( + x = mean(unlist(dataset[, input$indepvar])), y = mean(unlist(dataset[, input$outcome])), + name = "Center Point", marker = list(size = 20, color = "green", line = list(color = "yellow", width = 2)) + ) %>% + layout( + title = paste0( + "lm(", input$outcome, " ~ ", input$indepvar, + "), Cor(", input$indepvar, ",", input$outcome, ") = ", + round(cor(unlist(dataset[, input$indepvar]), unlist(dataset[, input$outcome])), 3) + ), + xaxis = list(title = input$indepvar), yaxis = list(title = input$outcome) + ) + }) + } +)