# 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.417318   95.619997
( 0.109165) ( 7.117342)

### $$\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)
 0.2438403
sd(DV)
 0.8985005

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 < 1e-4) { return(NA_real_) }
probs    <- diff(pnorm(brks, mean=param, sd=param))
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-Schätzer  0.3431997 1.1345665 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 < 1e-4) { return(NA_real_) } probs <- diff(pnorm(brks, mean=param, sd=param)) -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
 0.3541194 1.1531078