L[l] <- dista(route)
if(L[l] > L[l-1])         # worse tour
{
# energy
DL <- L[l] - L[l-1]
# transition probability
P_tr  <- exp(-DL/T)
u<- runif(1)
if(u <= P_tr)
{
# Metropolis: accepted
# generate new solution
r_sample<- sample(sL,replace = FALSE)
route<- c(s[1],r_sample,s[1])
}
else
{
# Metropolis: not accepted
L[l] <- L[l-1]
# generate new solution
r_sample<- sample(route0,replace = FALSE)
route<- c(s[1],r_sample,s[1])
}
} # ending worse tour
else
{
# better tour
best.route <- route
if (dista(best.route)<dista(very.best.route)) very.best.route <- best.route
# generate new solution
r_sample<- sample(sL,replace = FALSE)    # better tour
route<- c(s[1],r_sample,s[1])
if (dista(route)<dista(very.best.route)) very.best.route <- route
}
}    # ending tour length loop
minL[t]<- min(L)
}    # ending cooling loop
# Graphically show the route
# swap Long ang Lat
tmp <- cx
coords <- cx
coords[,1] <- tmp[,2]
coords[,2] <- tmp[,1]
xmin <- 38
xmax <- 46.5
ymin <- 7
ymax <- 17
# plot TSP
plot(seq(ymin,ymax,(ymax-ymin)/10),seq(xmin,xmax,(xmax-xmin)/10),type="n",xlab="Longitude",ylab="Latitude")
text(coords[s[1],2],coords[s[1],1],labels=labs[s[1]],col="black",cex=0.8)
text(coords[s[1],2],coords[s[1],1],labels="O",col="black",cex=2.5)
for (i in (2:20))
text(coords[s[i],2],coords[s[i],1],labels=labs[s[i]],col="red",cex=0.8)
# best route found by simulated annealing
for (i in 1:19)
arrows(coords[very.best.route[s[i]],2],coords[very.best.route[s[i]],1],coords[very.best.route[s[i+1]],2],
coords[very.best.route[s[i+1]],1],col="blue",angle=20,length=0.1)
arrows(coords[very.best.route[s[20]],2],coords[very.best.route[s[20]],1],coords[very.best.route[s[1]],2],
coords[very.best.route[s[1]],1],col="blue",angle=20,length=0.1)
# Show solution and route length
print(very.best.route)
print(dista(very.best.route))
# Code_8_4
# Read the province coordinates and mutual distances
setwd("C:/RPA/code/markov_chain_montecarlo")
city.coords <- read.table("data/provinces.txt",h=T)
cx <- as.matrix(city.coords[,3:4])
labs <- city.coords$Code
distances <- read.table("data/distances.txt",h=TRUE)
M_Cities <- data.matrix(distances)
# distance function
dista <- function(x) {
ret <- M_Cities[x[1],x[2]]
for (i in 2:20)
ret <- ret + M_Cities[x[i],x[i+1]]
as.numeric(ret)}
# Simulated Annealing (SA) code
r<-nrow(M_Cities) # number of raws of  M_Cities
# route starting from L'Aquila
s<-1:r
route <- c(s,1)
# route coming from SA
best.route <- route
# route corresponding to best minimum
very.best.route <- route
# Initial route length
L0 <- dista(route)
# Other parameters
# initial "temperature"
T0<-1000
# number of cooling steps
nT<-100
# number of tentative routes, for every T value
nL<-20
# slow-cooling parameter
a<-0.9
# route vector
L<-vector()
sL<-vector()
# best routes found during cooling
minL<-vector()
# Comment the following for an actual SA algorithm
set.seed(4321)
for (t in 1:nT)
{  # starting cooling loop
T<-a^t*T0
L[1]<-L0
r_sample<-sample(s[2:20],replace = FALSE)
sL<-r_sample
route<- c(s[1],r_sample,s[1])
for (l in 2:nL)
{  # starting tour length loop
route0<- route[2:20]
L[l] <- dista(route)
if(L[l] > L[l-1])         # worse tour
{
# energy
DL <- L[l] - L[l-1]
# transition probability
P_tr  <- exp(-DL/T)
u<- runif(1)
if(u <= P_tr)
{
# Metropolis: accepted
# generate new solution
r_sample<- sample(sL,replace = FALSE)
route<- c(s[1],r_sample,s[1])
}
else
{
# Metropolis: not accepted
L[l] <- L[l-1]
# generate new solution
r_sample<- sample(route0,replace = FALSE)
route<- c(s[1],r_sample,s[1])
}
} # ending worse tour
else
{
# better tour
best.route <- route
if (dista(best.route)<dista(very.best.route)) very.best.route <- best.route
# generate new solution
r_sample<- sample(sL,replace = FALSE)    # better tour
route<- c(s[1],r_sample,s[1])
if (dista(route)<dista(very.best.route)) very.best.route <- route
}
}    # ending tour length loop
minL[t]<- min(L)
}    # ending cooling loop
# Graphically show the route
# swap Long ang Lat
tmp <- cx
coords <- cx
coords[,1] <- tmp[,2]
coords[,2] <- tmp[,1]
xmin <- 38
xmax <- 46.5
ymin <- 7
ymax <- 17
# plot TSP
plot(seq(ymin,ymax,(ymax-ymin)/10),seq(xmin,xmax,(xmax-xmin)/10),type="n",xlab="Longitude",ylab="Latitude")
text(coords[s[1],2],coords[s[1],1],labels=labs[s[1]],col="black",cex=0.8)
text(coords[s[1],2],coords[s[1],1],labels="O",col="black",cex=2.5)
for (i in (2:20))
text(coords[s[i],2],coords[s[i],1],labels=labs[s[i]],col="red",cex=0.8)
# best route found by simulated annealing
for (i in 1:19)
arrows(coords[very.best.route[s[i]],2],coords[very.best.route[s[i]],1],coords[very.best.route[s[i+1]],2],
coords[very.best.route[s[i+1]],1],col="blue",angle=20,length=0.1)
arrows(coords[very.best.route[s[20]],2],coords[very.best.route[s[20]],1],coords[very.best.route[s[1]],2],
coords[very.best.route[s[1]],1],col="blue",angle=20,length=0.1)
# Show solution and route length
print(very.best.route)
print(dista(very.best.route))
setwd("C:/RPA/code/bayesian")
## Code_9_6.R
# Arrival time
# Total observation time
t<-120
# Arrival times for bus A
# Frequency of bus A
lambda.A<-1/10
set.seed(456)
# Number of arrivals of bus A
N.A<-rpois(1,lambda.A*t)
unifs<-runif(N.A,0,t)
arrivals.A<-sort(unifs)
plot(arrivals.A,rep(1,length(arrivals.A)),yaxt="n",ylab="",xlab="Bus arrival times",ylim=c(0,4),xlim=c(0,120),pch="A",cex=0.9)
# Arrival times for bus B
# Frequency of bus B
lambda.B<-1/15
set.seed(789)
# Number of arrivals of bus B
N.B<-rpois(1,lambda.B*t)
unifs<-runif(N.B,0,t)
arrivals.B<-sort(unifs)
points(arrivals.B,rep(2,length(arrivals.B)),pch="B",cex=0.9)
# Arrival times for bus C
# Frequency of bus C
lambda.C<-1/20
set.seed(123)
# Number of arrivals of bus C
N.C<-rpois(1,lambda.C*t)
unifs<-runif(N.C,0,t)
arrivals.C<-sort(unifs)
points(arrivals.C,rep(3,length(arrivals.C)),pch="C",cex=0.9)
## Code_9_6.R
# Arrival times
# Total observation time
t <- 120
#Arrival times for bus A
lambda<-1/10
set.seed(123)
N<-rpois(1,lambda*t)
unifs<-runif(N,0,t)
(arrivals.A<-sort(unifs))
#Arrival times for bus B
lambda<-1/15
set.seed(987)
N<-rpois(1,lambda*t)
unifs<-runif(N,0,t)
(arrivals.B<-sort(unifs))
#Arrival times for bus C
lambda<-1/20
set.seed(543)
N<-rpois(1,lambda*t)
unifs<-runif(N,0,t)
arrivals.C<-sort(unifs)
N1 <- length(arrivals.A)
T1 <- arrivals.A[N1]
N2 <- length(arrivals.B)
T2 <- arrivals.B[N2]
N3 <- length(arrivals.C)
T3 <- arrivals.C[N3]
Jeffreys (improper) Prior = 1/lambda
Posterior = Gamma(n,T)
yi = lambda_i
set.seed(111)
lam1 <- rgamma(50000,N1,T1)
set.seed(222)
lam2 <- rgamma(50000,N2,T2)
set.seed(333)
lam3 <- rgamma(50000,N3,T3)
hist(lam1,freq=FALSE,main="",xlab=expression(lambda[A]))
lines(seq(0,0.25,0.001),dgamma(seq(0,0.25,0.001),N1,T1))
hist(lam2,freq=FALSE,main="",xlab=expression(lambda[B]))
lines(seq(0,0.25,0.001),dgamma(seq(0,0.25,0.001),N2,T2))
hist(lam3,freq=FALSE,main="",xlab=expression(lambda[C]),ylim=c(0,15))
lines(seq(0,0.25,0.001),dgamma(seq(0,0.25,0.001),N3,T3))
c(N1/T1,mean(lam1))
c(N2/T2,mean(lam2))
c(N3/T3,mean(lam3))
p <- lam2/(lam1+lam2+lam3)
(h <- hist(p,freq=FALSE,main="",xlab=expression(p(lambda[B]/paste(Sigma,lambda[i])))))
mean(p)
setwd("C:/RPA/code/genetic_algorithm")
setwd("C:/RPA/code/genetic_algorithm")
## Code_10_4.R
# Read monthly runoff
setwd("C:/RPA/code/genetic_algorithm")
dat <- read.table("data/loire_runoff.txt",h=TRUE)
names(dat)
x <- ts(dat$DISCHRG)
# annual average
Y <- dat$YEAR
N <- length(dat$YEAR)
years <- Y[1]:Y[N]
n <- Y[N]-Y[1]+1
X <- rep(0,n)
for (i in 1:n){
X[i] <- sum(x[(12*(i-1)+1):(12*(i-1)+12)])
}
X[113] <- 0.5*(X[112]+X[114])
# Standardized annual runoff
X.m <- mean(X)
X.sd <- sd(X)
X <- ts((X-X.m)/X.sd,start=1863)
plot(X,xlab="Year",ylab=expression(paste("Annual runoff (m"^"3","/s)")),mgp=c(2,0.7,0))
# Count the number of observations
n <- length(X)
# Initialize error
e <- rep(1, n)
# Fitness function
fitness.fun <- function(crom){
g <- numeric
mu <- crom[1]
P <- round(crom[2])
Q <- round(crom[3])
if (P>Q)
for (i in 1:P) e[i] <- 0
else
e[1] <- 0
# Take into account a non-zero mean value
XX <- X-mu
fi <- c()
teta <- c()
if (P>0)
for (i in 1:P) fi = cbind(fi,crom[i+3])
if (Q>0)
for (j in 1:Q) teta = cbind(teta,crom[j+P+3])
if (P==0)
n_start <- 2
else
n_start <- P+1
for (t in (n_start : n)){
X_hat <- 0
if (P>0){
for (i in 1:P)
X_hat = X_hat+fi[i]*XX[t-i]
}
if (Q>0 & Q<=P){
for (i in 1:Q)
X_hat = X_hat+teta[i]*e[t-i]
}
e[t] <- (XX[t] - X_hat)
}
g <-  e%*%e
return(-g-(P+Q+1))
}
# Utility function for postprocessing
ARMA.params <- function(crom){
mu <- crom[1]
p <- round(crom[2])
q <- round(crom[3])
fi <- c()
teta <- c()
if (p>0)
for (i in 1:p) fi = cbind(fi,crom[i+3])
if (q>0)
for (j in 1:q) teta = cbind(teta,crom[j+p+3])
cat("mean",mu,"\n")
cat("AR",p,fi,"\n")
cat("MA",q,teta,"\n")
}
# Load GA library
library(GA)
# Initialize GA
pop.size<- 30
max.gen <- 300
xover.prob <- 0.7
mut.prob <- 0.1
# Compute
GA <- ga(type = "real-valued", fitness = fitness.fun,
lower = c(-0.01,0,0,-0.9,-0.9,-0.9,-0.9),
upper=c(0.01,2,2,0.9,0.9,0.9,0.9),
maxiter=max.gen,popSize = pop.size,elitism=4)
# Show solution
sol <- GA@solution[1,]
ARMA.params(sol)
# Plot diagnostics
plot(GA)
## Code_10_2.R
# Non linear fitting with GA
library(GA)
# Initialize
pop.size<- 100
chrome.len <- 60
max.gen <- 100
xover.prob <- 0.7
mut.prob <- 0.1
# Read data
setwd("C:/RPA/code/genetic_algorithm")
dat <- read.table("data/exp_data.txt")
X <- dat$V1
Y <- dat$V2
# Define fitness function
f <- function(x)  {
chrom.a <- x[1:20]
a <- a1+(a2-a1)*binary2decimal(chrom.a)/dmax
chrom.b <- x[21:40]
b <- b1+(b2-b1)*binary2decimal(chrom.b)/dmax
chrom.c <- x[41:60]
c <- c1+(c2-c1)*binary2decimal(chrom.c)/dmax
N <- length(X)
ret <- 0
for (i in 1:N)
ret <- ret + (Y[i] - a*exp(-b*X[i])+c)^2
1/ret
}
# Exec GA
a1 <- 0
a2 <- 1.0
b1 <- 0
b2 <- 1.0
c1 <- 0
c2 <- 1.0
dmax = binary2decimal(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1))
GA <- ga(type = "binary", fitness = f, popSize = pop.size, nBits = chrome.len, maxiter = max.gen, pmutation = mut.prob, pcrossover = xover.prob,keepBest = TRUE)
result <- unname(GA@solution)
a <- a1+(a2-a1)*binary2decimal(result[1:20])/dmax
b <- b1+(b2-b1)*binary2decimal(result[21:40])/dmax
c <- c1+(c2-c1)*binary2decimal(result[41:60])/dmax
param <- c(a,b,c)
# Output
max(GA@fitness)
min(GA@fitness)
mean(GA@fitness)
param
# Plot
# Exponential function
fexp <- function(x,p)
{
a <- p[1]
b <- p[2]
c <- p[3]
a*exp(-b*x)+c
}
plot(X,Y)
curve(fexp(x,param),from=0,to=10,add=TRUE,col="red",lty=2)
# Diagnostics
plot(GA@summary[,1],t="l",xlab="Generation",ylab="Fitness")
## Code_10_5.R
library(GA)
set.seed(12345)
# N cities: 1 2 .. N
N <- 15
# city coordinates
coords <- matrix(nrow=N,ncol=2)
test <- TRUE
while (test) {
coords[,1] <- trunc(10*runif(N))
coords[,2] <- trunc(10*runif(N))
# avoid coincident cities
test <- anyDuplicated(coords)
}
distances <- matrix(rep( 0, len=N*N), nrow = N)
for (i in 1:N){
for (j in 1:N){
distances[i,j] <- sqrt((coords[i,1]-coords[j,1])^2+(coords[i,2]-coords[j,2])^2)
}
}
# fitness function
inverse_distance <- function(x) {
N <- length(x)
dist.total <- distances[1,x[1]]
for (i in 1:(N-1))
dist.total <- dist.total + distances[x[i],x[i+1]]
dist.total <- dist.total + distances[x[N],1]
1/dist.total
}
TSP.ga <- ga(type = "permutation", fitness = inverse_distance, lower = 2, upper = N,
elitism = 1, maxiter = 500, popSize = 100)
# show GA parameters
summary(TSP.ga)
route.best <- unname(TSP.ga@solution)
# choose the first one
route.best <- c(1, route.best[1,])
# plot TSP
plot(0:10,0:10,type="n",xlab="",ylab="")
text(coords[1,1],coords[1,2],labels=1,col="black")
text(coords[1,1],coords[1,2],labels="O",col="black",cex=2)
for (i in (2:N))
text(coords[i,1],coords[i,2],labels=i,col="red")
# best route find by GA
for (i in 1:(N-1))
arrows(coords[route.best[i],1],coords[route.best[i],2],coords[route.best[i+1],1],
coords[route.best[i+1],2],col="blue",angle=20,length=0.1)
arrows(coords[route.best[N],1],coords[route.best[N],2],coords[route.best[1],1],
coords[route.best[1],2],col="blue",angle=20,length=0.1)
# plot algorithm convergence
plot(TSP.ga)
## Code_10_5.R
library(GA)
set.seed(12345)
# N cities: 1 2 .. N
N <- 15
# city coordinates
coords <- matrix(nrow=N,ncol=2)
test <- TRUE
while (test) {
coords[,1] <- trunc(10*runif(N))
coords[,2] <- trunc(10*runif(N))
# avoid coincident cities
test <- anyDuplicated(coords)
}
distances <- matrix(rep( 0, len=N*N), nrow = N)
for (i in 1:N){
for (j in 1:N){
distances[i,j] <- sqrt((coords[i,1]-coords[j,1])^2+(coords[i,2]-coords[j,2])^2)
}
}
# fitness function
inverse_distance <- function(x) {
N <- length(x)
dist.total <- distances[1,x[1]]
for (i in 1:(N-1))
dist.total <- dist.total + distances[x[i],x[i+1]]
dist.total <- dist.total + distances[x[N],1]
1/dist.total
}
TSP.ga <- ga(type = "permutation", fitness = inverse_distance, lower = 2, upper = N,
elitism = 1, maxiter = 500, popSize = 100)
# show GA parameters
summary(TSP.ga)
route.best <- unname(TSP.ga@solution)
# choose the first one
route.best <- c(1, route.best[1,])
# plot TSP
plot(0:10,0:10,type="n",xlab="",ylab="")
text(coords[1,1],coords[1,2],labels=1,col="black")
text(coords[1,1],coords[1,2],labels="O",col="black",cex=2)
for (i in (2:N))
text(coords[i,1],coords[i,2],labels=i,col="red")
# best route find by GA
for (i in 1:(N-1))
arrows(coords[route.best[i],1],coords[route.best[i],2],coords[route.best[i+1],1],
coords[route.best[i+1],2],col="blue",angle=20,length=0.1)
arrows(coords[route.best[N],1],coords[route.best[N],2],coords[route.best[1],1],
coords[route.best[1],2],col="blue",angle=20,length=0.1)
# plot algorithm convergence
plot(TSP.ga)
