# Numerical Optimization

## Use cases

### Maximum-likelihood estimation

Weibull distribution

library(MASS)
X <- rweibull(100, shape=1.5, scale=100)
fitdistr(X, densfun="weibull", start=list(shape=1, scale=50))
     shape        scale
1.4418009   90.8713787
( 0.1142759) ( 6.5847708)

### $$\chi^{2}$$ test of normality

The $$\chi^{2}$$ normality test that compares observed and expected category frequencies after partitioning the value range into discrete bins (categories) does not work correctly when the sample mean and standard deviation are taken as estimates $$\hat{\mu}$$ and $$\hat{\sigma}$$. The test statistic then does not actually have a $$\chi^{2}$$ distribution.

DV <- rnorm(50, mean=0, sd=1)
mean(DV)
[1] -0.09821994
sd(DV)
[1] 0.9775147

Instead, it is necessary to estimate $$\hat{\mu}$$ and $$\hat{\sigma}$$ taking into account the chosen categories. One possibility to do this is by finding the minimum $$\chi^{2}$$ estimates, i.e., those estimates that minimize the observed $$X^{2}$$ statistic.

nCls <- 4  # number of classes

## category limits -> cut points for equal-probability intervals
limits <- qnorm(seq(from=1/nCls, to=(nCls-1)/nCls, length.out=nCls-1),
mean=0, sd=1)

## determine observed category frequencies
breaks   <- c(-Inf, limits, Inf)
DVcut    <- cut(DV, breaks=breaks)
observed <- table(DVcut)

## function that calculates the empirical X^2 test statistic
## needs to be minimized -> objective function
minFunMin <- function(param, brks, obs) {
if(param[2] < 1e-4) { return(NA_real_) }
probs    <- diff(pnorm(brks, mean=param[1], sd=param[2]))
expected <- sum(obs) * probs
sum((obs-expected)^2 / expected)
}

Optimization = minimization of objective function

resMinChisq <- optim(c(mean(DV), sd(DV)), minFunMin,
brks=breaks, obs=observed, gr=NULL, method="BFGS")
resMinChisq$par # min-chi^2-estimate [1] -0.05177565 0.99517343 The other possibility is a grouped maximum likelihood estimate - where the likelihood comes from the multinomial distribution. ## objective function ## negative multinomial log likelihood minFunGML <- function(param, brks, obs) { if(param[2] < 1e-4) { return(NA_real_) } probs <- diff(pnorm(brks, mean=param[1], sd=param[2])) -dmultinom(obs, size=sum(obs), prob=probs, log=TRUE) } ## optimize -> minimize objective function resGrML <- optim(c(mean(DV), sd(DV)), minFunGML, brks=breaks, obs=observed, gr=NULL, method="BFGS") resGrML$par
[1] -0.05208171  0.99805581