# Ouloumuc June 2012 Tutorial 01 ### Bernoulli(m,n,p) random context #subsets are identified to binary vectors #Bernoulli context BC = function(rows,col,pr) { BC = matrix(rbinom(rows*col,1,pr),rows,col) } # Galois connection (f,g) #subset of items g = function(itemset, context) { y = rep(1,m) #if itemset is not empty if (1 %in% itemset) { ones = which(itemset ==1) for (j in ones) { y = y*context[,j] } } g = y } #subset of objects f = function(objset, context) { y = rep(1,n) #if set of objects is not empty if (1 %in% objset) { ones = which(objset == 1) for (i in ones) { y = y*context[i,] } } f = y } #Closure functions h = function(objset, context) { h = g(f(objset,context),context)} k = function(itemset, context) { k = f(g(itemset,context),context)} #Closed sets kclosed = function(itemset,context) { kclosed = (min(itemset == k(itemset,context)))} hclosed = function(objset,context) { hclosed = (min(objset == h(objset,context)))} #Simulating sim Random Bernoulli(p) Context #itemset empirical probability i_eprob = function(sim,A) { s = 0 for (i in 1:sim) { D = BC(m,n,p) s = s + kclosed(A,D) } i_eprob = s/sim } #object set empirical probability o_eprob = function(sim,A) { s = 0 for (i in 1:sim) { D = BC(m,n,p) s = s + hclosed(A,D) } o_eprob = s/sim } ##non binary vectors c(1,2,...) #convert a vector c(1,2,...) in a binary one v2b = function(vect,K) { y = rep(0,K) l = length(vect) #if vect is not the empty vector if (l > 0) { for (i in vect) { y[i] = 1 } } v2b = y } #convert a binary into c(1,2,...) b2v = function(binv,K) { y = c() #if vect is not the empty vector if (1 %in% binv) { y = which(binv == 1) } b2v = y } #Galois connection in non binary Gf = function(vector, context) { Gf = b2v(f(v2b(vector,m),context),n) } Gg = function(vector, context) { Gg = b2v(f(v2b(vector,n),context),m) } Gh = function(vector, context) { Gh = b2v(h(v2b(vector,m),context),m) } Gk = function(vector, context) { Gk = b2v(k(v2b(vector,n),context),n) } ### Exact probability that an itemset be closed pr_ci = function(size, m0,n0,p0){ s = 0 for (k in 0:m0) { s = s+ choose(m0,k)*p0^(k*size)*(1- p^(size))^(m0-k)*(1 - p^k)^(n0 - size) } pr_ci = s } #Example #probability of 1 p = 0.4 # size of the context m = 15 n = 7 #number of simulations s = 10000 #Is itemset {1,4} is k-closed ? Sometimes yes, Sometimes no. D = BC(m,n,p) r = Gk(c(1,4,5),D) D r #itemset IS = c(1,2,5,7) I = v2b(IS,n) # empirical and exact probability ep = i_eprob(s,I) pr = pr_ci(length(IS),m,n,p) print(paste("Empirical probability:",ep)) print(paste("Exact probability:",pr)) #Object set OS = c(2,6,10) O = v2b(OS,m) ep = o_eprob(s,O) ep ### simulating a Markov chain state = c('A','C','G','T') l = length(state) p0 = rgamma(l,1,1) p0 =p0/sum(p0) index = 1:l X=c() X[1] = sample(index,1, prob = p0) MM = matrix(rgamma(l^2,1,1),l,l) s = apply(MM,1,'sum') for (i in 1:l) { MM[i,]=MM[i,]/s[i] } sim = 15 for (i in 1:sim) { X[i+1] = sample(index,1, prob = MM[X[i],]) } X = state[X] X DELETE 1