forked from Jpomz/detecting-spectra-differences
/
MLE_tidy.R
49 lines (45 loc) · 1.47 KB
/
MLE_tidy.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
MLE_tidy <- function(df, rsp_var){
# define variables
x <- df[[rsp_var]]
xmin = min(x)
xmax = max(x)
log.x = log(x)
sum.log.x = sum(log.x)
# initial starting point for parameter estimate
PL.bMLE = 1/(log(min(x)) - sum.log.x/length(x)) - 1
# non-linear minimization
PLB.minLL = nlm(negLL.PLB,
p = PL.bMLE, x = x, n = length(x),
xmin = xmin, xmax = xmax,
sumlogx = sum.log.x)
# estimate for b
PLB.bMLE = PLB.minLL$estimate
# minimum estimate of b
PLB.minNegLL = PLB.minLL$minimum
## 95% CI calculation
bvec = seq(PLB.bMLE - 0.5, PLB.bMLE + 0.5, 1e-05)
PLB.LLvals = vector(length = length(bvec))
for (i in 1:length(bvec)) {
PLB.LLvals[i] = negLL.PLB(bvec[i],
x = x,
n = length(x),
xmin = xmin,
xmax = xmax,
sumlogx = sum.log.x)
}
critVal = PLB.minNegLL + qchisq(0.95, 1)/2
bIn95 = bvec[PLB.LLvals < critVal]
# confidence interval
PLB.MLE.bConf = c(min(bIn95), max(bIn95))
if (PLB.MLE.bConf[1] == min(bvec) |
PLB.MLE.bConf[2] == max(bvec)) {
dev.new()
plot(bvec, PLB.LLvals)
abline(h = critVal, col = "red")
stop("Need to make bvec larger - see R window")
}
# return b estimate and min/max 95% CI
return(data.frame(b = PLB.bMLE,
minCI = min(bIn95),
maxCI = max(bIn95)))
}