File: reliability.R

package info (click to toggle)
r-cran-markovchain 0.8.5-4-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 2,060 kB
  • sloc: cpp: 2,854; sh: 13; makefile: 2
file content (152 lines) | stat: -rw-r--r-- 4,830 bytes parent folder | download | duplicates (3)
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
#load required libraries
library(parallel)
require(MCMCpack)
require(markovchain)
dimensions2Test<-2:32
numSim=10000 

#helper function to create a random stochastic matrix

createMatrix<-function(matr_size) {
  out<-matrix(0, nrow=matr_size, ncol = matr_size)
  for (i in 1:matr_size) {
    priors.dirichlet<-runif(n=matr_size) 
    myStochasticRow<-rdirichlet(n=1,alpha=priors.dirichlet)
    out[i,]<-myStochasticRow
  }
  return(out)
}


createSparseMatrix<-function(matr_size, sparsity=0.75){
  out <- matrix(0, nrow=matr_size, ncol = matr_size)
  nonzeroitems<-ceiling(matr_size*(1-sparsity))
  for (i in 1:matr_size) {
    priors.dirichlet<-runif(n=nonzeroitems)
    myStochasticRow<-rdirichlet(n=1,alpha=priors.dirichlet)
    columnsPositions<-sample(x = 1:matr_size,size = nonzeroitems,replace = FALSE)
    out[i,columnsPositions]<-myStochasticRow
  }
  return(out)
}


#first test: function to simulate the inversion of a matrix of size num
checkInversion<-function(i,num){
  #simulate the priors
  myStochasticMatrix<-createMatrix(matr_size = num)
  #this code returns FALSE -> 0 if error in inversion 1 otherwise
  out<-tryCatch(steadyStates(as(myStochasticMatrix, "markovchain")),
                error=function(c) return(FALSE)
  )
  if(class(out)=="logical") return(0) else return(1)
}

checkSparseMInversion<-function(i,num){
  #simulate the priors
  myStochasticMatrix<-createSparseMatrix(matr_size = num)
  #this code returns FALSE -> 0 if error in inversion 1 otherwise
  out<-tryCatch(steadyStates(as(myStochasticMatrix, "markovchain")),
                error=function(c) return(FALSE)
  )
  if(class(out)=="logical") return(0) else return(1)
}




#performing the simulation

successRate<-numeric(length(dimensions2Test))

#using parallel backend
no_cores <- detectCores() - 1 
cl <- makeCluster(no_cores)
clusterExport(cl, "checkInversion")
clusterExport(cl, "createMatrix")
clusterEvalQ(cl, library(markovchain))
clusterEvalQ(cl, library(MCMCpack))
k=1
for (dimension in dimensions2Test){
  simulations<-parSapply(cl=cl,1:numSim,FUN=checkInversion,num=dimension)
  successRate[k]<-mean(simulations)
  k=k+1
}
stopCluster(cl)
#summarising first test:
#plot(x=dimensions2Test,y=successRate,type="l",xlab="matrix size",ylab="success rate",main="Steady state computation success rate")
#abline(h=0.5,col="red")
#text(x=dimensions2Test,y=successRate,labels=round(successRate,digits=2),col="darkred",cex=0.7)
#dev.off()

dimensions2Test = 2^seq(from=3, to=8)
successRate<-numeric(length(dimensions2Test))

#using parallel backend
no_cores <- detectCores() - 1 
cl <- makeCluster(no_cores)
clusterExport(cl, "checkSparseMInversion")
clusterExport(cl, "createSparseMatrix")
clusterEvalQ(cl, library(markovchain))
clusterEvalQ(cl, library(MCMCpack))
k=1
for (dimension in dimensions2Test){
  simulations<-parSapply(cl=cl,1:numSim,FUN=checkSparseMInversion,num=dimension)
  successRate[k]<-mean(simulations)
  k=k+1
}
stopCluster(cl)
# 
# 
# plot(x=dimensions2Test,y=successRate,type="l",xlab="matrix size",ylab="success rate",main="Steady state computation success rate ??? sparse matrices")
# abline(h=0.5,col="red")
# text(x=dimensions2Test,y=successRate,labels=round(successRate,digits=2),col="darkred",cex=0.7)



#second test: simulating exponentiation

checkExponentiation<-function(i,num){
  #simulate the priors
  myStochasticMatrix<-createMatrix(matr_size = num)
  #this code returns FALSE -> 0 if error in inversion 1 otherwise
  out<-tryCatch((as(myStochasticMatrix, "markovchain"))^2,
                error=function(c) return(FALSE)
  )
  if(class(out)=="logical") return(0) else return(1)
}


#performing the simulation
successRate2<-numeric(length(dimensions2Test))


#using parallel backend
no_cores <- detectCores() - 1 
cl <- makeCluster(no_cores)
clusterExport(cl, "checkExponentiation")
clusterExport(cl, "createMatrix")
clusterEvalQ(cl, library(markovchain))
clusterEvalQ(cl, library(MCMCpack))
k=1
for (dimension in dimensions2Test){
  simulations<-parSapply(cl=cl,1:numSim,FUN=checkExponentiation,num=dimension)
  successRate2[k]<-mean(simulations)
  k=k+1
}
stopCluster(cl)
#summarising first test:


#par(mfrow=c(1,2))

# plot(x=dimensions2Test,y=successRate,type="l",xlab="matrix size",ylab="success rate",main="Steady state computation success rate")
# abline(h=0.5,col="red")
# text(x=dimensions2Test,y=successRate,labels=round(successRate,digits=2),col="darkred",cex=0.7)



# plot(x=dimensions2Test,y=successRate2,type="l",xlab="matrix sixe",ylab="success rate",main="Exponentiation computation success rate")
# abline(h=0.5,col="red")
# text(x=dimensions2Test,y=successRate2,labels=round(successRate2,digits=2),col="darkred",cex=0.7)