Wednesday, September 29, 2004
Hi, David. I was just wondering, are you going to post the code that you went over in class last night? Or do you have some other samples of code that runs in R? I find that looking at examples of code is a good way to learn what it can do, and how to write it. Thanks!!
Here's the R code:
Here's the R code:
# read the data
myData <- read.table("c:\\download\\fatalityOrdered.txt")
x <- myData[,1]
y <- myData[,2]
# constants
N <- dim(myData)[1]
MAXWINDOWSIZE <- 250
MINWINDOWSIZE <- 5
NUMITER <- 10
# initialization
pValue <- 0;
smallestLambdaArray <- rep(0,NUMITER);
count <- 1 #first iteration of the algorithm is to find actual clusters, all other permutations are for comparison
compareprob <- rep(0, NUMITER)
while(count <= NUMITER) {
cat("iteration: ",count,"\n");
smallestLambda <- 100000
windowSize <- MINWINDOWSIZE
while(windowSize < MAXWINDOWSIZE){
windowStart <- 1
while (windowStart <= (N - windowSize + 1)) {
window <- y[windowStart:(windowStart+windowSize-1)]
positives <- sum(window)
negatives <- windowSize - positives
outWindowSize <- N - windowSize
outPositives <- sum(y) - positives
outNegatives <- sum(1-y) - negatives
if ((positives > 0) && (negatives > 0) && ((positives/windowSize)>(sum(y)/N))) {
windowLambda <- -1 * (log((positives/windowSize)^positives) + log((negatives/windowSize)^negatives) + log((outPositives/outWindowSize)^outPositives) + log((outNegatives/outWindowSize)^outNegatives))
if (windowLambda < smallestLambda) {
smallestLambda <- windowLambda
# cat(x[windowStart],"-",x[windowStart+windowSize-1],"\t",windowLambda,"\t",windowSize,"\t",positives,"\n");
if(count == 1) {
windowSizeOfCluster <- windowSize
clusterWindowStart <- windowStart
}
}
}
windowStart <- windowStart + 1
}
windowSize <- windowSize + 1
}
y <- sample(y, N, replace = FALSE)
smallestLambdaArray[count] <- smallestLambda
if (smallestLambdaArray[count] <= smallestLambdaArray[1]) {
pValue <- pValue + 1;
}
count <- count + 1
}
cat("Size: ", windowSizeOfCluster,"\n")
cat("Window: ",x[clusterWindowStart],"-",x[clusterWindowStart+windowSizeOfCluster-1],"\n")
cat("pValue: ",pValue/NUMITER,"\n")
