1### R code from vignette source 'setpartitions.Rnw' 2 3################################################### 4### code chunk number 1: load.library 5################################################### 6 7 8 9################################################### 10### code chunk number 2: setpartitions.Rnw:82-83 11################################################### 12require(partitions) 13 14 15################################################### 16### code chunk number 3: alleles 17################################################### 18setparts(c(2,2,1)) 19 20 21################################################### 22### code chunk number 4: use_achims_idea 23################################################### 24a <- c(1,2,3,2,1) 25split(seq_along(a),a) 26 27 28################################################### 29### code chunk number 5: setpartitions.Rnw:362-365 30################################################### 31options(width=63) 32m <- 9 33n <- 3 34 35 36################################################### 37### code chunk number 6: setpartsrestrictedparts 38################################################### 39jj <- setparts(restrictedparts(m,n,include.zero=FALSE)) 40summary(jj) 41 42 43################################################### 44### code chunk number 7: definepenaltyfunctionF 45################################################### 46tau <- 1:9 47slowest <- function(x) max(tapply(tau, x, sum)) 48 49 50################################################### 51### code chunk number 8: show_the_best 52################################################### 53time.taken <- apply(jj,2,slowest) 54minimal.time <- sum(tau)/n 55jj[,time.taken == minimal.time] 56 57 58################################################### 59### code chunk number 9: define.E.matrix 60################################################### 61E <- matrix(c(0,0,0,1,1,0,0,0,1,1,1,1,1,0,0), 5,3) 62dimnames(E) <- list( 63 evidence=paste("E",1:5,sep=""), 64 crime=paste("C",1:3,sep="") 65 ) 66 67 68################################################### 69### code chunk number 10: print.E 70################################################### 71E 72 73 74################################################### 75### code chunk number 11: set.evidence.matrix 76################################################### 77a <- 78structure( 79 c(2, 2, 4, 4, 3, 1, 2, 1, 1, 2, 2, 1, 4, 2, 1, 2, 1, 80 4, 2, 4, 2, 1, 4, 4, 1, 2, 1, 4, 2, 1, 1, 2, 1, 1, 2), 81 .Dim = c(5L, 7L), 82 .Dimnames = list( 83 evidence=paste("E",1:5,sep=""), 84 crime=paste("C",1:7,sep="") 85 ) 86) 87 88 89################################################### 90### code chunk number 12: show.a 91################################################### 92a 93 94 95################################################### 96### code chunk number 13: calc.likelihoods 97################################################### 98genbeta <- function(x){gamma(sum(x))/prod(gamma(x))} 99 100lp <- function(evidence, alpha, partition){ 101 r <- length(unique(partition)) 102 evidence <- as.factor(evidence) 103 levels(evidence) <- seq_along(alpha) 104 out <- 1 105 for(k in unique(partition)){ #k is a perp 106 thisperp <- partition==k 107 evidence.thisperp <- evidence[thisperp] 108 no.of.crimes.thisperp <- sum(thisperp) 109 out <- out/genbeta(alpha+table(evidence.thisperp)) 110 } 111 return(out*genbeta(alpha)^r) 112} 113 114sp <- setparts(7) 115l1 <- apply(sp,2, function(x){lp(evidence=a[1,],alpha=rep(1,2),partition=x)}) 116l2 <- apply(sp,2, function(x){lp(evidence=a[2,],alpha=rep(1,2),partition=x)}) 117l3 <- apply(sp,2, function(x){lp(evidence=a[3,],alpha=rep(1,4),partition=x)}) 118l4 <- apply(sp,2, function(x){lp(evidence=a[4,],alpha=rep(1,4),partition=x)}) 119l5 <- apply(sp,2, function(x){lp(evidence=a[5,],alpha=rep(1,4),partition=x)}) 120 121likelihood <- l1*l2*l3*l4*l5 122likelihood <- likelihood/max(likelihood) 123support <- log(likelihood) 124 125 126################################################### 127### code chunk number 14: plotlikelihoods 128################################################### 129par(mfcol=c(1,2)) 130plot(likelihood,ylab="likelihood",xlab="hypothesis") 131abline(h = exp(-2)) 132plot(support,ylab="support",xlab="hypothesis") 133abline(h = -2) 134 135 136################################################### 137### code chunk number 15: show.maximum.likelihood 138################################################### 139sp <- setparts(7) 140 141 142################################################### 143### code chunk number 16: setpartitions.Rnw:673-674 144################################################### 145sp[,which.max(support)] 146 147 148################################################### 149### code chunk number 17: likelihood.gt.2 150################################################### 151index <- order(support,decreasing=TRUE) 152sp <- sp[,index] 153support <- support[index] 154dimnames(sp) <- list( 155 crime = paste("C",1:7, sep=""), 156 partition = paste("H",1:ncol(sp),sep="") 157 ) 158 159 160################################################### 161### code chunk number 18: setpartitions.Rnw:695-696 162################################################### 163 sp[, support > -2] 164 165 166################################################### 167### code chunk number 19: show.likelihoods 168################################################### 169support[support > -2] 170 171 172