1
"estimateTruncatedNormalParameters" <-
2
function(numericalSample, testLevel = 0.975){
4
# Checking for NAs in sample
5
numericalSample <- numericalSample[ ! is.na(numericalSample)]
7
a <- min(numericalSample)
8
b <- max(numericalSample)
9
# First guess for the optimization step. Should find a better solution.
10
mean <- mean(numericalSample)
11
var <- var(numericalSample)
16
phiA <- dnorm(a, mean, sd)
17
phiB <- dnorm(b, mean, sd)
18
PhiA <- pnorm(a, mean, sd)
19
PhiB <- pnorm(b, mean, sd)
20
denom <- (PhiB - PhiA - (mean0 - a) * phiA - (b - mean0) * phiB)
21
var <- var0 * (PhiB - PhiA) / denom
22
mean <- mean0 + var * (phiB - phiA) / (PhiB - PhiA)
25
# For reasons of data types and calls, we need to use an auxiliary function. However,
26
# this is just an intermediary step, we still use fitdistr on the Beta Pdf.
27
AuxiliaryTruncatedNormal <- function(numericalSample,...){
28
# Call to fitdistr (MASS package) .
29
estimation <- fitdistr(numericalSample, computeTruncatedNormalPdf,
30
start = list(mu = mean, sigma = sqrt(var)), ...)
34
try(estimation <- AuxiliaryTruncatedNormal(numericalSample, a = (a - .Machine$double.eps), b = (b + .Machine$double.eps)))
35
if (is.na(estimation))
37
# If the first method of estimation fails, we try the L-BFGS-B method.
38
estimation <- AuxiliaryTruncatedNormal(numericalSample, a = (a - .Machine$double.eps), b = (b + .Machine$double.eps), method = "L-BFGS-B", lower = c(-Inf,0))
41
# Checking whether loglikelihood is supported in the current R version.
42
if (length(estimation$loglik) == 0) {
43
warning("Warning: the current R version does not support the loglikelihood computation")
44
estimation$loglik = "Warning: the current R version does not support the loglikelihood computation"
47
# Computation of the t-value used in the CIs
49
tValue <- qt(testLevel, length(numericalSample)-1)
51
confidenceIntervalmu <- c(as.numeric(estimation$estimate[1]) -
52
tValue * as.numeric(estimation$sd[1]),
53
as.numeric(estimation$estimate[1]) +
54
tValue * as.numeric(estimation$sd[1]))
56
confidenceIntervalsigma <- c(as.numeric(estimation$estimate[2]) -
57
tValue * as.numeric(estimation$sd[2]),
58
as.numeric(estimation$estimate[2]) +
59
tValue * as.numeric(estimation$sd[2]))
61
estimatedTruncatedNormal <- list(distribution = "TruncatedNormal",
62
mu = as.numeric(estimation$estimate[1]),
63
sigma = as.numeric(estimation$estimate[2]),
66
confidenceIntervalmu = confidenceIntervalmu,
67
confidenceIntervalsigma = confidenceIntervalsigma,
68
logLikelihood = estimation$loglik)
added
removed
debian/rotRPackage/R/estimateTruncatedNormalParameters.R
