shiny app to give Edgeworth approximation to the density based on input raw cumulants.

shiny_edgeworth.R

library(mpoly)

# blinnikov and moessner's ADVANCE function
advance <- function(kms) {
	current <- kms$current
	mold <- kms$mold

	n <- length(current)
	ords <- 1:(length(current))
	stopifnot(n == sum(current * ords))
	sumcur <- n
	m <- 1
	is.done <- FALSE
	while (!is.done) {
		sumcur <- sumcur - current[m]*m + (m+1)
		current[m] <- 0
		current[m+1] <- current[m+1] + 1
		m <- m+1
		is.done <- (sumcur <= n) || (m > mold)
	}
	mold <- max(mold,m)
	current[1] <- n - sumcur
	retv <- list(current=current,mold=mold)
	return(retv)
}

dapx_edgeworth <- function(raw.cumulants) {#FOLDUP
	order.max <- length(raw.cumulants)

	mu <- raw.cumulants[1]
	sigma <- sqrt(raw.cumulants[2])

	retval <- mpoly(list(c(eta=0,coef=1)))

	# compute via equation (43) of Blinnikov and Moessner
	if (order.max > 2) {
		Sn <- raw.cumulants[3:order.max] / (raw.cumulants[2] ^ (2:(order.max-1)))

		for (s in c(1:(order.max-2))) {
			nexterm <- mpoly(list(c(eta=0,coef=0)))
			kms <- list(mold=1,current=c(s,rep(0,s-1)))
			while (kms$mold <= s) {
				r <- sum(kms$current)
				coefs <- ((Sn[1:s] / factorial(3:(s+2))) ^ kms$current) / factorial(kms$current)
				coef <- prod(coefs)
				nexterm <- nexterm + coef * plug(mpoly::hermite(degree=s+2*r+1, normalized=FALSE),'x','eta')
				kms <- advance(kms)
			}
			retval <- retval + (sigma^s) * nexterm
		}
	}
	# but what is eta, really?
	eta <- (1/sigma) * mpoly(list(c(x=1,coef=1),c(x=0,coef=-mu)))
	# plug in, and adjust back from standardized
	retval <- (1/sigma) * plug(retval, "eta", eta)
	return(retval)
}#UNFOLD


library(shiny)

server <- function(input,output) {
	output$polynomial_rep <- renderText({
		raw.cumulants <- c(input$mean,input$stdev^2,input$cumul3,input$cumul4,input$cumul5)
		mu <- raw.cumulants[1]
		sigma <- sqrt(raw.cumulants[2])
		eta <- (1/sigma) * mpoly(list(c(x=1,coef=1),c(x=0,coef=-mu)))

		mypoly <- dapx_edgeworth(raw.cumulants)
		retval <- paste0("phi(",print(eta,stars=TRUE),") * (",print(mypoly,stars=TRUE),")")
		retval
	})
}

ui <- fluidPage(
	sidebarLayout(
    sidebarPanel(
      numericInput("mean","mean:",value=0),
      numericInput("stdev","stdev:",value=1,min=0),
      numericInput("cumul3","raw 3rd cumulant:",value=1),
      numericInput("cumul4","raw 4th cumulant:",value=0),
      numericInput("cumul5","raw 5th cumulant:",value=0)
    ),
    mainPanel(textOutput("polynomial_rep"))
  )
)
shinyApp(ui=ui,server=server)