Np <- length(dd)
# Needed matrices
Hj1 <- matrix(0,N,Np)
Hj2 <- matrix(0,N,Np)
G = matrix(0,2,N)
hj1 <- matrix(0,Np)
hj2 <- matrix(0,Np)
hh <- matrix(0,Np)
PP <- matrix(0,Np)
ss <- matrix(0,Np)
# Decompose in the orthogonal base sin/cos
k <- 1
m <- 2
for (d in dd)
{
w <- 4*pi*d*BW/(N*c)
C <- N/2 + sin(N*w)/(2*sin(w))
S <- N/2 - sin(N*w)/(2*sin(w))
# Gjk matrix
for (i in 1:N)
{
l <- -T + 2*T*(i-1)/(N-1)
G[1,i] <- cos(w*l)
G[2,i] <- sin(w*l)
}
# orthogonal functions
Hj1[,k] <- G[1,]/sqrt(C)
Hj2[,k] <- G[2,]/sqrt(S)
# project data on the orthogonal base
h1 <- sum(si*Hj1[,k]);
h2 <- sum(si*Hj2[,k]);
hj1[k] <- h1;
hj2[k] <- h2;
h2m = 0.5*(h1^2+h2^2)
hh[k] = h2m;
# Periodogram: equation (9.30)
P <- (1-m*h2m/(N*d2m))^((m-N)/2)
PP[k] <- P
# Periodogram: equation (9.36)
sigma2 <- N/(N-m-2)*(d2m - m*h2m/N)
#P <- exp(m*h2m/sigma2)
ss[k] <- sigma2
#PP[k] <- P
k <- k+1
}
HHMAX <- max(hh)
SSMIN <- min(ss)
# Posterior probability: equation (9.30)
p <- PP/max(PP)
# Plots
plot(dd,p,t="l",xlab="Distance (m)",ylab="Posterior Probability",mgp=c(2,0.8,0))
plot(dd,abs(hh),t="l",xlab="Distance (m)",ylab=expression(paste("<",h^2,">")),mgp=c(2,0.8,0))
# Signal to Noise Ratio
snr <- (2*HHMAX/N)/SSMIN
## Code_9_7.R
setwd("C:/RPA/code/bayesian")
# speed of EM signal in air
c <- 2.998e8
# read data: multiple measurements
M <- readLines("data/multiple_rain.txt")
temp <- list()
vec <- c(list())
for (i in 1:length(M))
{
one.line <- M[i]
strsplit(one.line,"SIGNAL,") -> tmp
header <- tmp[[1]][[1]]
M.textdata <- paste(header,"SIGNAL,")
dat <- tmp[[1]][[2]]
dat <- as.numeric(unlist(strsplit(dat, ",")))
M.data <- dat
temp$data <- M.data
temp$textdata <- M.textdata
vec[[i]] <- temp
}
# vec[[1]][[1]], vec[[2]][[1]] etc contengono i dati
# vec[[1]][[2]], vec[[2]][[2]] etc contengono le intestazioni
choice <- 40
meas <- vec[[as.integer(choice)]][[1]];
# Elaboration
f1 <- 76e9
f2 <- 77e9
N <- length(meas)
deltaf <- (f2-f1)/(N-1)
f <- f1 + (0:(N-1))*deltaf
BW <- f2-f1
# normalized signal
meas <- meas-mean(meas)
si <- meas/max(meas)
# Bayesian spectrum analysis
d2m <- mean(si^2)
T <- round((N-1)/2)
T1 <- -T
T2 <- T-1
# choose distance interval
dd <- seq(2,9,0.01)
Np <- length(dd)
# Needed matrices
Hj1 <- matrix(0,N,Np)
Hj2 <- matrix(0,N,Np)
G = matrix(0,2,N)
hj1 <- matrix(0,Np)
hj2 <- matrix(0,Np)
hh <- matrix(0,Np)
PP <- matrix(0,Np)
ss <- matrix(0,Np)
# Decompose in the orthogonal base sin/cos
k <- 1
m <- 2
for (d in dd)
{
w <- 4*pi*d*BW/(N*c)
C <- N/2 + sin(N*w)/(2*sin(w))
S <- N/2 - sin(N*w)/(2*sin(w))
# Gjk matrix
for (i in 1:N)
{
l <- -T + 2*T*(i-1)/(N-1)
G[1,i] <- cos(w*l)
G[2,i] <- sin(w*l)
}
# orthogonal functions
Hj1[,k] <- G[1,]/sqrt(C)
Hj2[,k] <- G[2,]/sqrt(S)
# project data on the orthogonal base
h1 <- sum(si*Hj1[,k]);
h2 <- sum(si*Hj2[,k]);
hj1[k] <- h1;
hj2[k] <- h2;
h2m = 0.5*(h1^2+h2^2)
hh[k] = h2m;
# Periodogram: equation (9.30)
P <- (1-m*h2m/(N*d2m))^((m-N)/2)
PP[k] <- P
# Periodogram: equation (9.36)
sigma2 <- N/(N-m-2)*(d2m - m*h2m/N)
#P <- exp(m*h2m/sigma2)
ss[k] <- sigma2
#PP[k] <- P
k <- k+1
}
HHMAX <- max(hh)
SSMIN <- min(ss)
# Posterior probability: equation (9.30)
p <- PP/max(PP)
# Plot periodogram
plot(dd,abs(hh),t="l",xlab="Distance (m)",ylab=expression(paste("<",h^2,">")),mgp=c(2,0.8,0))
# Signal to Noise Ratio
snr <- (2*HHMAX/N)/SSMIN
## Code_9_8.R
setwd("C:\RPA\bayesian\code")
# speed of EM signal in air
c <- 2.998e8
# read data: multiple measurements
M <- readLines("data/multiple_rain.txt")
temp <- list()
vec <- c(list())
for (i in 1:length(M))
{
one.line <- M[i]
strsplit(one.line,"SIGNAL,") -> tmp
header <- tmp[[1]][[1]]
M.textdata <- paste(header,"SIGNAL,")
dat <- tmp[[1]][[2]]
dat <- as.numeric(unlist(strsplit(dat, ",")))
M.data <- dat
temp$data <- M.data
temp$textdata <- M.textdata
vec[[i]] <- temp
}
# vec[[1]][[1]], vec[[2]][[1]] etc contengono i dati
# vec[[1]][[2]], vec[[2]][[2]] etc contengono le intestazioni
# average measurements
Nstart <- 40
Nave <- 5
meas = 0;
for (i in 1:Nave)
meas <- meas + vec[[Nstart+i-1]][[1]]
meas <- meas/Nave
# Elaboration
f1 <- 76e9
f2 <- 77e9
N <- length(meas)
deltaf <- (f2-f1)/(N-1)
f <- f1 + (0:(N-1))*deltaf
BW <- f2-f1
# normalized signal
meas <- meas-mean(meas)
si <- meas/max(meas)
# Bayesian spectrum analysis
d2m <- mean(si^2)
T <- round((N-1)/2)
T1 <- -T
T2 <- T-1
# choose distance interval
dd <- seq(2,9,0.01)
Np <- length(dd)
# Needed matrices
Hj1 <- matrix(0,N,Np)
Hj2 <- matrix(0,N,Np)
G = matrix(0,2,N)
hj1 <- matrix(0,Np)
hj2 <- matrix(0,Np)
hh <- matrix(0,Np)
PP <- matrix(0,Np)
ss <- matrix(0,Np)
# Decompose in the orthogonal base sin/cos
k <- 1
m <- 2
for (d in dd)
{
w <- 4*pi*d*BW/(N*c)
C <- N/2 + sin(N*w)/(2*sin(w))
S <- N/2 - sin(N*w)/(2*sin(w))
# Gjk matrix
for (i in 1:N)
{
l <- -T + 2*T*(i-1)/(N-1)
G[1,i] <- cos(w*l)
G[2,i] <- sin(w*l)
}
# orthogonal functions
Hj1[,k] <- G[1,]/sqrt(C)
Hj2[,k] <- G[2,]/sqrt(S)
# project data on the orthogonal base
h1 <- sum(si*Hj1[,k]);
h2 <- sum(si*Hj2[,k]);
hj1[k] <- h1;
hj2[k] <- h2;
h2m = 0.5*(h1^2+h2^2)
hh[k] = h2m;
# Periodogram: equation (9.30)
P <- (1-m*h2m/(N*d2m))^((m-N)/2)
PP[k] <- P
# Periodogram: equation (9.36)
sigma2 <- N/(N-m-2)*(d2m - m*h2m/N)
P <- exp(m*h2m/sigma2)
ss[k] <- sigma2
PP[k] <- P
k <- k+1
}
HHMAX <- max(hh)
SSMIN <- min(ss)
# Posterior probability: equation (9.30)
p <- PP/max(PP)
# Plot periodogram
plot(dd,abs(hh),t="l",xlab="Distance (m)",ylab=expression(paste("<",h^2,">")),mgp=c(2,0.8,0))
# Signal to Noise Ratio
snr <- (2*HHMAX/N)/SSMIN
## Code_9_8.R
setwd("C:\RPA\bayesian\code")
# speed of EM signal in air
c <- 2.998e8
# read data: multiple measurements
M <- readLines("data/multiple_rain.txt")
setwd("C:/RPA/code/bayesian")
## Code_9_8.R
setwd("C:/RPA/code/bayesian")
## Code_9_8.R
setwd("C:/RPA/code/bayesian")
# speed of EM signal in air
c <- 2.998e8
# read data: multiple measurements
M <- readLines("data/multiple_rain.txt")
temp <- list()
vec <- c(list())
for (i in 1:length(M))
{
one.line <- M[i]
strsplit(one.line,"SIGNAL,") -> tmp
header <- tmp[[1]][[1]]
M.textdata <- paste(header,"SIGNAL,")
dat <- tmp[[1]][[2]]
dat <- as.numeric(unlist(strsplit(dat, ",")))
M.data <- dat
temp$data <- M.data
temp$textdata <- M.textdata
vec[[i]] <- temp
}
# vec[[1]][[1]], vec[[2]][[1]] etc contengono i dati
# vec[[1]][[2]], vec[[2]][[2]] etc contengono le intestazioni
# average measurements
Nstart <- 40
Nave <- 5
meas = 0;
for (i in 1:Nave)
meas <- meas + vec[[Nstart+i-1]][[1]]
meas <- meas/Nave
# Elaboration
f1 <- 76e9
f2 <- 77e9
N <- length(meas)
deltaf <- (f2-f1)/(N-1)
f <- f1 + (0:(N-1))*deltaf
BW <- f2-f1
# normalized signal
meas <- meas-mean(meas)
si <- meas/max(meas)
# Bayesian spectrum analysis
d2m <- mean(si^2)
T <- round((N-1)/2)
T1 <- -T
T2 <- T-1
# choose distance interval
dd <- seq(2,9,0.01)
Np <- length(dd)
# Needed matrices
Hj1 <- matrix(0,N,Np)
Hj2 <- matrix(0,N,Np)
G = matrix(0,2,N)
hj1 <- matrix(0,Np)
hj2 <- matrix(0,Np)
hh <- matrix(0,Np)
PP <- matrix(0,Np)
ss <- matrix(0,Np)
# Decompose in the orthogonal base sin/cos
k <- 1
m <- 2
for (d in dd)
{
w <- 4*pi*d*BW/(N*c)
C <- N/2 + sin(N*w)/(2*sin(w))
S <- N/2 - sin(N*w)/(2*sin(w))
# Gjk matrix
for (i in 1:N)
{
l <- -T + 2*T*(i-1)/(N-1)
G[1,i] <- cos(w*l)
G[2,i] <- sin(w*l)
}
# orthogonal functions
Hj1[,k] <- G[1,]/sqrt(C)
Hj2[,k] <- G[2,]/sqrt(S)
# project data on the orthogonal base
h1 <- sum(si*Hj1[,k]);
h2 <- sum(si*Hj2[,k]);
hj1[k] <- h1;
hj2[k] <- h2;
h2m = 0.5*(h1^2+h2^2)
hh[k] = h2m;
# Periodogram: equation (9.30)
P <- (1-m*h2m/(N*d2m))^((m-N)/2)
PP[k] <- P
# Periodogram: equation (9.36)
sigma2 <- N/(N-m-2)*(d2m - m*h2m/N)
P <- exp(m*h2m/sigma2)
ss[k] <- sigma2
PP[k] <- P
k <- k+1
}
HHMAX <- max(hh)
SSMIN <- min(ss)
# Posterior probability: equation (9.30)
p <- PP/max(PP)
# Plot periodogram
plot(dd,abs(hh),t="l",xlab="Distance (m)",ylab=expression(paste("<",h^2,">")),mgp=c(2,0.8,0))
# Signal to Noise Ratio
snr <- (2*HHMAX/N)/SSMIN
## Code_9_6.R
# Arrival times
# 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_11.R
## Where is the light ?
# distance from sea
Y0 <- 1
# distance along the coast
X0 <- 2
# possible values of X (positions of photo-detectors)
dx <- 0.05
X <- seq(-5,5,dx)
Nx <- length(X)
######################################################################
# number of measurements
# to change to observe the effect on the distribution a posteriori
N <- 10
######################################################################
# tetak <- runif(k,-pi/2,pi/2)
# instead of taking theta between -pi/2 and pi/2
# theta is selected such to determine "possible" x
# included between -x_max and x_max
x.max <- 50
# Here the angle is obtained using the arctan function, where x.max in an angle
tetak.max <- atan(x.max)
tetak <- runif(N,-tetak.max,tetak.max)
tetak
# compute the positions of the detectors activated by flash light
xk <- X0+Y0*tan(tetak)
# What is the distribution of the positions ?
hist(xk,main="")
L <- rep(0,Nx)
for (i in 1:Nx){
lk <- log(Y0^2+(xk-X[i])^2)
L[i] <- sum(lk)
}
hist(lk)
# posteriori probability
post <- exp(-L)
plot(X,dx*post/sum(post),t="l",ylab="p(X|x,Y0)")
abline(v=2,col="blue",lty=2)
abline(v=mean(xk),col="red")
## Code_Append_C_1.R
# Use rjags library
library(rjags)
# Export JAGS model
model <- textConnection("model {
Y ~ dbin(0.5,10)
Pr <- step(3.1 - Y)
}")
# Init
jags.inits <- NULL
# Data
jags.data <- list()
# Perform Bayesian analysis using JAGS
model <- jags.model(model, data=jags.data, inits=jags.inits, n.chains=1)
update(model, n.iter=10000)
# Parameters
jags.params <- c("Y","Pr") # parameters to be monitored
samples <- coda.samples(model, jags.params,10000)
# Analysis of the simulation results
plot(samples)
summary(samples)
# extract P and Y
out <- do.call(rbind.data.frame, samples)
names(out)
# Probability
summary(out$Pr)
sum(out$Pr==0)/length(out$Pr)
sum(out$Pr==1)/length(out$Pr)
# Flipping results
summary(out$Y)
hist(out$Y,breaks=c(0,1,2,3,4,5,6,7,8,9,10),main="",freq=FALSE,xlab="Number of heads")
## Code_Append_C_1.R
# Use rjags library
library(rjags)
install.packages("jagsUI")
## Code_Append_C_1.R
# Use rjags library
library(rjags)
install.packages("rjags")
## Code_Append_C_1.R
# Use rjags library
library(rjags)
## Code_Append_C_1.R
# Use rjags library
library(rjags)
# Export JAGS model
model <- textConnection("model {
Y ~ dbin(0.5,10)
Pr <- step(3.1 - Y)
}")
# Init
jags.inits <- NULL
# Data
jags.data <- list()
# Perform Bayesian analysis using JAGS
model <- jags.model(model, data=jags.data, inits=jags.inits, n.chains=1)
update(model, n.iter=10000)
# Parameters
jags.params <- c("Y","Pr") # parameters to be monitored
samples <- coda.samples(model, jags.params,10000)
# Analysis of the simulation results
plot(samples)
summary(samples)
# extract P and Y
out <- do.call(rbind.data.frame, samples)
names(out)
# Probability
summary(out$Pr)
sum(out$Pr==0)/length(out$Pr)
sum(out$Pr==1)/length(out$Pr)
# Flipping results
summary(out$Y)
hist(out$Y,breaks=c(0,1,2,3,4,5,6,7,8,9,10),main="",freq=FALSE,xlab="Number of heads")
## Code_Append_C_1.R
# Use rjags library
library(rjags)
install.packages("rjags")
## Code_Append_C_1.R
# Use rjags library
library(rjags)
## Code_Append_C_1.R
# Use rjags library
library(rjags)
