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