library(expm)
## Loading required package: Matrix
##
## Attaching package: 'expm'
## The following object is masked from 'package:Matrix':
##
## expm
library(markovchain)
## Package: markovchain
## Version: 0.6.9.14
## Date: 2019-01-20
## BugReport: http://github.com/spedygiorgio/markovchain/issues
library(diagram)
## Loading required package: shape
library(pracma)
##
## Attaching package: 'pracma'
## The following objects are masked from 'package:expm':
##
## expm, logm, sqrtm
## The following objects are masked from 'package:Matrix':
##
## expm, lu, tril, triu
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
#
#Shembuj te punuar Leksione Procese Rastid he Simulim ne R
# Pergatiti: PhD. Eralda Gjika (Dhamo)
# Departamenti i Matematikes se Aplikuar
# Fakulteti i Shkencave te Natyres
# Universiteti i Tiranes
# Email: eralda.dhamo@fshn.edu.al
# LinkedIn:https://www.linkedin.com/in/eralda-dhamo-gjika-71879128/?originalSubdomain=al
#
#
# USHTRIMI SHEMBULL Leksione Eralda Gjika
# Jepet Zinxhiri i Markovit me hapesrie gjendjesh E dhe matrice te Probabilitetet te kalimit Q
E<- c("SPSS","Matlab","R")# hapesira e gjendjeve
Q <- matrix(c(.65,.3,.05,.1,0.75,.15,0,.1,.9),nrow=3, byrow=TRUE)
row.names(Q) <- E; colnames(Q) <- E
Q # Matrica e probabiliteteve tekalimit me nje hap
## SPSS Matlab R
## SPSS 0.65 0.30 0.05
## Matlab 0.10 0.75 0.15
## R 0.00 0.10 0.90
# Ndertimi i grafit te probabiliteteve te kalimit me nje hap
plotmat(Q,pos = c(1,2),
lwd = 1, box.lwd = 2,
cex.txt = 0.8,
box.size = 0.1,
box.type = "circle",
box.prop = 0.5,
box.col = "light yellow",
arr.length=.1,
arr.width=.1,
self.cex = .4,
self.shifty = -.01,
self.shiftx = .13
)
Q2 <- Q %^% 2 # Matrica e probabiliteteve te kalimit me dy hapa
round(Q2,3) # rrumbullakos deri ne tre shifra pas presjes dhjetore vlerat e matrices Q2
## SPSS Matlab R
## SPSS 0.452 0.425 0.122
## Matlab 0.140 0.608 0.252
## R 0.010 0.165 0.825
f0<- c(0.5,0.2,0.3)# Deklarojme vektorin e probabiliteteve fillestare
f2=round(f0 %*% Q2,3) # situata probabilitare pas 2 momenteve te kohes
f2
## SPSS Matlab R
## [1,] 0.257 0.384 0.359
Q3 <- Q %^% 3 # llogaritja e matrices se probabiliteteve te kalimit me tre hapa
Q3 # Afishon matricen Q3
## SPSS Matlab R
## SPSS 0.336625 0.466750 0.196625
## Matlab 0.151750 0.522875 0.325375
## R 0.023000 0.209250 0.767750
round(Q3,3)# Matrica e rendit 3 me elemnte te rrumbullakosur deri ne tre shifra pas presjes dhjetore
## SPSS Matlab R
## SPSS 0.337 0.467 0.197
## Matlab 0.152 0.523 0.325
## R 0.023 0.209 0.768
f3=round(f2 %*% Q,3)# llogaritja e f3 nepermjet f2
f3
## SPSS Matlab R
## [1,] 0.205 0.401 0.394
round(f0 %*% Q3,3)# llogaritja e f3 nepermjet f0
## SPSS Matlab R
## [1,] 0.206 0.401 0.394
Q5 <- Q %^% 5# Matrica e rendit 5
round(Q5,3)
## SPSS Matlab R
## SPSS 0.220 0.459 0.321
## Matlab 0.145 0.436 0.419
## R 0.047 0.264 0.689
# shiko me poshte funksionin matrix.power(Q,fuqia) na ndimon te njehsojme
# matricene probabiliteteve t ekalimit me n-hapa
matrix.power <- function(mat, n)
{
if (n == 1) return(mat)
result <- diag(1, ncol(mat))
while (n > 0) {
if (n %% 2 != 0) {
result <- result %*% mat
n <- n - 1
}
mat <- mat %*% mat
n <- n / 2
}
return(result)
}
# si mund te perfitohet matrica Q5 nga matricat Q2 dhe Q3 per lehtesi veprimesh
Q2=round(matrix.power(Q,2),5)
Q3=round(matrix.power(Q,3),5)
Q2%*%Q3
## SPSS Matlab R
## [1,] 0.21963180 0.4590615 0.3213064
## [2,] 0.14512242 0.4358302 0.4190520
## [3,] 0.04737995 0.2635740 0.6890477
Q5
## SPSS Matlab R
## SPSS 0.2196341 0.4590594 0.3213066
## Matlab 0.1451231 0.4358272 0.4190497
## R 0.0473800 0.2635731 0.6890469
Q50=round(matrix.power(Q,50),3)
Q50
## SPSS Matlab R
## [1,] 0.098 0.341 0.561
## [2,] 0.098 0.341 0.561
## [3,] 0.098 0.341 0.561
# Si do te ishte situata probabiliteteve te kalimit grafikisht pas 5 momenteve
plotmat(round(Q5,3),pos = c(1,2),
lwd = 1, box.lwd = 2,
cex.txt = 0.8,
box.size = 0.1,
box.type = "circle",
box.prop = 0.5,
box.col = "light yellow",
arr.length=.1,
arr.width=.1,
self.cex = .4,
self.shifty = -.01,
self.shiftx = .13
)
#Krahasjme dy grafet e probabiliteteve te kalimit
par(mfrow=c(1,2))
# grafi per Q
plotmat(Q,pos = c(1,2),
lwd = 1, box.lwd = 2,
cex.txt = 0.8,
box.size = 0.1,
box.type = "circle",
box.prop = 0.5,
box.col = "light yellow",
arr.length=.1,
arr.width=.1,
self.cex = .4,
self.shifty = -.01,
self.shiftx = .13
)
# grafi per Q5
plotmat(round(Q5,3),pos = c(1,2),
lwd = 1, box.lwd = 2,
cex.txt = 0.8,
box.size = 0.1,
box.type = "circle",
box.prop = 0.5,
box.col = "light yellow",
arr.length=.1,
arr.width=.1,
self.cex = .4,
self.shifty = -.01,
self.shiftx = .13
)
#
Q10 <- Q %^% 10# Matrica e rendit 10
round(Q10,3)
## SPSS Matlab R
## SPSS 0.130 0.386 0.484
## Matlab 0.115 0.367 0.518
## R 0.081 0.318 0.600
Q50 <- Q %^% 50# Matrica e rendit 50
round(Q50,3)
## SPSS Matlab R
## SPSS 0.098 0.341 0.561
## Matlab 0.098 0.341 0.561
## R 0.098 0.341 0.561
#
#
####################
#Grafiku i konvergjences me 10 hapa
Hapi_Q <- function(f0, Q, n) {
f_n <- matrix(NA, n+1, length(f0))
f_n[1,] <- f0
for (i in seq_len(n))
f_n[i+1,] <- f0 <- f0 %*% Q
f_n
}
### GRAFI I KONVERGJENCES
y1=Hapi_Q(c(1,0,0),Q,10)
y2=Hapi_Q(c(0,1,0),Q,10)
y3=Hapi_Q(c(0,0,1),Q,10)
matplot(0:10, y1, type="l", lty=1, xlab="Hapi", ylab="y", las=1,main="Grafiku i konvergjences me 10 hapa")
matlines(0:10, y2, type="l", lty=2)
matlines(0:10, y3, type="l", lty=3)
