-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfitGeneralMMUsingCV.R
159 lines (108 loc) · 6.19 KB
/
fitGeneralMMUsingCV.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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
## 6.22.10
# Use cross validation to choose the proper # of mixture components for general mixture models
# Mixture model is fitted to data via EM -- see "emForGeneralMixtureModel.R"
# Randomly split data into k folds
splitDataIntoKFolds = function(dataVector, k) {
n = length(dataVector)
randomOrdering=sample(x=seq(from=1,to=n,by=1), size=n, prob=rep(1/n,n))
numPointsPerFold = floor(n/k)
foldAssignmentVector = vector(mode="numeric", length=n)
for (i in 1:n) {
positionInRandomOrdering = which(randomOrdering == i)
foldAssignmentVector[i] = ceiling(positionInRandomOrdering/numPointsPerFold)
}
# Allocate any remaining points if n is not evenly divisible by k
extraPoints = which(foldAssignmentVector == k+1)
numExtraPoints = length(extraPoints)
if (numExtraPoints > 0) {
for (i in 1:numExtraPoints) {
foldAssignmentVector[extraPoints[i]] = i
}
}
return(foldAssignmentVector)
}
fitMMToTrainingDataAndReturnValidationLikelihood = function(trainingDataVector, validationDataVector, numMixtureComponents, distrType, initType) {
#source("/Users/radlab/Desktop/ksauer/Desktop/scads/experiments/client/performance/logparsing/src/main/R/emFitAndSample.R")
# likelihoodValues = vector(mode="numeric", length=numRdmInits)
# for (i in 1:numRdmInits) {
# emParams = fitEM(observationsVector=trainingDataVector, k=numMixtureComponents, distrType="gaussian", modelName="V")
# likelihoodValues[i] = computeGMMLogLikelihood(dataVector=validationDataVector, parameters=emParams)
# }
# if (aggregationFunction == "mean") {
# return(mean(likelihoodValues))
# } else if (aggregationFunction == "median") {
# return(median(likelihoodValues))
# } else if (aggregationFunction == "max") {
# return(max(likelihoodValues))
# } else {
# print("Unsupported aggregation function.")
# }
source("emForGeneralMixtureModel.R")
# fit MM to training data
emFit = em(trainingDataVector, numMixtureComponents, distrType, initType)
# get loglik of validation data under MM
validationLoglik = computeLoglik(validationDataVector, emFit$params, distrType)
# return validation loglik
return(validationLoglik)
}
findAvgLogLikelihoodViaCVFixedNumMixtureComponents = function(dataVector, numFolds, numMixtureComponents, distrType, initType, destinationPath="") {
# Partition the data into folds
partitionVector = splitDataIntoKFolds(dataVector=dataVector, k=numFolds)
# Get log likelihood for each fold
# i => validation set
# all other folds => training set
logLikelihood = vector(mode="numeric", length=numFolds)
for (i in 1:numFolds) {
print(paste("Using fold ", i, " for validation data", sep=""))
validationData = dataVector[which(partitionVector == i)]
trainingData = dataVector[which(partitionVector != i)]
#logLikelihood[i] = fitGMMToTrainingDataAndReturnValidationLikelihood(trainingDataVector=trainingData, validationDataVector=validationData, numMixtureComponents=numMixtureComponents, numRdmInits=numRdmInits, aggregationFunction="median")
logLikelihood[i] = fitMMToTrainingDataAndReturnValidationLikelihood(trainingData, validationData, numMixtureComponents, distrType, initType)
}
if (destinationPath != "") {
save(logLikelihood, file=paste(destinationPath, "/cvLogLikelihood-", numMixtureComponents, "mixtureComponents.RData", sep=""))
}
# Avg the log likelihood over all folds
return(mean(logLikelihood))
}
kFoldCVToChooseNumMixtureComponentsForMM = function(dataVector, distrType, initType, destinationPath, numFolds=10, mixtureComponentsRange=seq(from=2,by=1,to=10)) {
dir.create(destinationPath)
avgLogLikelihood = vector(mode="numeric", length=length(mixtureComponentsRange))
for (i in 1:length(mixtureComponentsRange)) {
print(paste("# mixture components =", mixtureComponentsRange[i]))
#avgLogLikelihood[i] = findAvgLogLikelihoodViaCVFixedNumMixtureComponents(dataVector=dataVector, numFolds=numFolds, numMixtureComponents=mixtureComponentsRange[i], numRdmInits=numRdmInits, destinationPath=destinationPath)
avgLogLikelihood[i] = findAvgLogLikelihoodViaCVFixedNumMixtureComponents(dataVector, numFolds, mixtureComponentsRange[i], distrType, initType, destinationPath)
}
return(mixtureComponentsRange[which.max(avgLogLikelihood)])
}
plotCVLogLikelihoodVsNumMixtureComponents = function(path, mixtureComponentsRange=seq(from=2,by=1,to=10)) {
pdf(file=paste(path, "/cvLogLikelihoodVsNumMixtureComponents.pdf", sep=""))
par(mar=c(5,5,4,2)+0.1)
# Computing stats for each run
minPerNumMixtureComponents = vector(mode="numeric", length=length(mixtureComponentsRange))
maxPerNumMixtureComponents = vector(mode="numeric", length=length(mixtureComponentsRange))
avgPerNumMixtureComponents = vector(mode="numeric", length=length(mixtureComponentsRange))
medianPerNumMixtureComponents = vector(mode="numeric", length=length(mixtureComponentsRange))
for (i in 1:length(mixtureComponentsRange)) {
load(paste(path, "/cvLogLikelihood-", mixtureComponentsRange[i], "mixtureComponents.RData", sep="")) # => logLikelihood
minPerNumMixtureComponents[i] = min(logLikelihood)
maxPerNumMixtureComponents[i] = max(logLikelihood)
avgPerNumMixtureComponents[i] = mean(logLikelihood)
medianPerNumMixtureComponents[i] = median(logLikelihood)
}
# Assumes loglikelihood values are < 0
ymin = 1.05*min(minPerNumMixtureComponents)
ymax = 0.95*max(maxPerNumMixtureComponents)
# Assumes loglikelihood values are > 0
#ymin = 0.95*min(minPerNumMixtureComponents)
#ymax = 1.05*max(maxPerNumMixtureComponents)
plot(x=0,y=0,xlim=c(min(mixtureComponentsRange),max(mixtureComponentsRange)),ylim=c(ymin,ymax), col=0, xlab="# Mixture Components", ylab="Log Likelihood", main="Log Likelihood vs. # Mixture Components")
for (i in 1:length(mixtureComponentsRange)) {
load(paste(path, "/cvLogLikelihood-", mixtureComponentsRange[i], "mixtureComponents.RData", sep="")) # => logLikelihood
points(x=rep(mixtureComponentsRange[i], length(logLikelihood)), y=logLikelihood, col="blue", lw=2)
}
lines(x=mixtureComponentsRange, y=avgPerNumMixtureComponents, lw=2, col="red")
lines(mixtureComponentsRange, medianPerNumMixtureComponents, lw=2, col="green")
legend("bottomright", legend=c("log likelihood per fold", "average", "median"), col=c("blue", "red", "green"), pch=c(1,-1,-1), lty=c(0,1,1), lwd=c(2,2,2))
dev.off()
}