.Rhistory

# p is used to generate snow load
p <- runif(length(id_snow))
# Final generation
g_s[id_snow] <- (B + 1.0 / A * (-log(-NS * log(1.0 - p_e + p_e * p))))^2*rooffactor(100,0.68,0.22)/
((Wind_velocity_50)^2*rooffactor(100,0.68,0.22))
g_s
mean(g_s)
D_wind <- cumsum(T_wind,g_s)
g_s <-
T_wind <- rexp(N, 1/0.03835)
D_wind <- cumsum(T_wind,g_s)
D_wind <- stepfun(cumsum(T_wind),g_s)
length(T_wind)
length(g_s)
N <- 100
NS <- 24
G_cov <- 0.17
G_mean <- (1+3.050*G_cov)/(1+2.592*G_cov)
A <- 1.282/(G_cov*G_mean)
B <- G_mean - (0.577/A)
# p_0 is the probability of no-snow in one year (10 segment)
p0 <- exp(-exp(A*B))
# p_e is the probability of snow in a segment
p_e <- 1.00 - exp(-1.00/NS*(exp(A*B)))
# parameters for distribution W_50
Wind_velocity_50 <- B + 1.0 / A * (-log(-NS * log(1.0 - p_e + p_e * 49/50)))
# generate N random number to decide whether it snows in one segment
r_n <- runif(N+1)
# generated snow laod
g_s <- rep(0,N+1)
# id that it snows in that segment. with r_n < p_e
# it does not snow in that winter segment if r_n > p_e
id_snow <- which (r_n <= p_e)
# p is used to generate snow load
p <- runif(length(id_snow))
# Final generation
g_s[id_snow] <- (B + 1.0 / A * (-log(-NS * log(1.0 - p_e + p_e * p))))^2*rooffactor(100,0.68,0.22)/
((Wind_velocity_50)^2*rooffactor(100,0.68,0.22))
g_s
T_wind <- rexp(N, 1/0.03835)
length(T_wind)
length(g_s)
D_wind <- stepfun(cumsum(T_wind),g_s)
t <- seq(from = 0, to = max(cumsum(T_wind)),
by = 0.01)
plot(t, D_wind(t),type = 'l',main = 'wind_load')
N <- 100
NS <- 12
G_cov <- 0.17
G_mean <- (1+3.050*G_cov)/(1+2.592*G_cov)
A <- 1.282/(G_cov*G_mean)
B <- G_mean - (0.577/A)
# p_0 is the probability of no-snow in one year (10 segment)
p0 <- exp(-exp(A*B))
# p_e is the probability of snow in a segment
p_e <- 1.00 - exp(-1.00/NS*(exp(A*B)))
# parameters for distribution W_50
Wind_velocity_50 <- B + 1.0 / A * (-log(-NS * log(1.0 - p_e + p_e * 49/50)))
# generate N random number to decide whether it snows in one segment
r_n <- runif(N+1)
# generated snow laod
g_s <- rep(0,N+1)
# id that it snows in that segment. with r_n < p_e
# it does not snow in that winter segment if r_n > p_e
id_snow <- which (r_n <= p_e)
# p is used to generate snow load
p <- runif(length(id_snow))
# Final generation
g_s[id_snow] <- (B + 1.0 / A * (-log(-NS * log(1.0 - p_e + p_e * p))))^2*rooffactor(100,0.68,0.22)/
((Wind_velocity_50)^2*rooffactor(100,0.68,0.22))
g_s
mean(g_s)
rooffactor(100,0.68,0.22)/rooffactor(100,0.68,0.22)
T_wind <- rexp(N, 1/0.03835)
length(T_wind)
length(g_s)
D_wind <- stepfun(cumsum(T_wind),g_s)
t <- seq(from = 0, to = max(cumsum(T_wind)),
by = 0.01)
plot(t, D_wind(t),type = 'l',main = 'wind_load')
mean(g_s)
NS <- 12
G_cov <- 0.108
G_mean <- (1+3.050*G_cov)/(1+2.592*G_cov)
A <- 1.282/(G_cov*G_mean)
B <- G_mean - (0.577/A)
# p_0 is the probability of no-snow in one year (10 segment)
p0 <- exp(-exp(A*B))
# p_e is the probability of snow in a segment
p_e <- 1.00 - exp(-1.00/NS*(exp(A*B)))
# parameters for distribution W_50
Wind_velocity_50 <- B + 1.0 / A * (-log(-NS * log(1.0 - p_e + p_e * 49/50)))
# generate N random number to decide whether it snows in one segment
r_n <- runif(N+1)
# generated snow laod
g_s <- rep(0,N+1)
# id that it snows in that segment. with r_n < p_e
# it does not snow in that winter segment if r_n > p_e
id_snow <- which (r_n <= p_e)
# p is used to generate snow load
p <- runif(length(id_snow))
# Final generation
g_s[id_snow] <- (B + 1.0 / A * (-log(-NS * log(1.0 - p_e + p_e * p))))^2*rooffactor(100,0.68,0.22)/
((Wind_velocity_50)^2*rooffactor(100,0.68,0.22))
g_s
mean(g_s)
N <- 100
NS <- 12
G_cov <- 0.150
G_mean <- (1+3.050*G_cov)/(1+2.592*G_cov)
A <- 1.282/(G_cov*G_mean)
B <- G_mean - (0.577/A)
# p_0 is the probability of no-snow in one year (10 segment)
p0 <- exp(-exp(A*B))
# p_e is the probability of snow in a segment
p_e <- 1.00 - exp(-1.00/NS*(exp(A*B)))
# parameters for distribution W_50
Wind_velocity_50 <- B + 1.0 / A * (-log(-NS * log(1.0 - p_e + p_e * 49/50)))
# generate N random number to decide whether it snows in one segment
r_n <- runif(N+1)
# generated snow laod
g_s <- rep(0,N+1)
# id that it snows in that segment. with r_n < p_e
# it does not snow in that winter segment if r_n > p_e
id_snow <- which (r_n <= p_e)
# p is used to generate snow load
p <- runif(length(id_snow))
# Final generation
g_s[id_snow] <- (B + 1.0 / A * (-log(-NS * log(1.0 - p_e + p_e * p))))^2*rooffactor(100,0.68,0.22)/
((Wind_velocity_50)^2*rooffactor(100,0.68,0.22))
g_s
mean(g_s)
N <- 100
NS <- 4
G_cov <- 0.150
G_mean <- (1+3.050*G_cov)/(1+2.592*G_cov)
A <- 1.282/(G_cov*G_mean)
B <- G_mean - (0.577/A)
# p_0 is the probability of no-snow in one year (10 segment)
p0 <- exp(-exp(A*B))
# p_e is the probability of snow in a segment
p_e <- 1.00 - exp(-1.00/NS*(exp(A*B)))
# parameters for distribution W_50
Wind_velocity_50 <- B + 1.0 / A * (-log(-NS * log(1.0 - p_e + p_e * 49/50)))
# generate N random number to decide whether it snows in one segment
r_n <- runif(N+1)
# generated snow laod
g_s <- rep(0,N+1)
# id that it snows in that segment. with r_n < p_e
# it does not snow in that winter segment if r_n > p_e
id_snow <- which (r_n <= p_e)
# p is used to generate snow load
p <- runif(length(id_snow))
# Final generation
g_s[id_snow] <- (B + 1.0 / A * (-log(-NS * log(1.0 - p_e + p_e * p))))^2*rooffactor(100,0.68,0.22)/
((Wind_velocity_50)^2*rooffactor(100,0.68,0.22))
g_s
mean(g_s)
N <- 100
NS <- 4
G_cov <- 0.150
G_mean <- (1+3.050*G_cov)/(1+2.592*G_cov)
A <- 1.282/(G_cov*G_mean)
B <- G_mean - (0.577/A)
# p_0 is the probability of no-snow in one year (10 segment)
p0 <- exp(-exp(A*B))
# p_e is the probability of snow in a segment
p_e <- 1.00 - exp(-1.00/NS*(exp(A*B)))
# parameters for distribution W_50
Wind_velocity_50 <- B + 1.0 / A * (-log(-NS * log(1.0 - p_e + p_e * 49/50)))
# generate N random number to decide whether it snows in one segment
r_n <- runif(N+1)
# generated snow laod
g_s <- rep(0,N+1)
# id that it snows in that segment. with r_n < p_e
# it does not snow in that winter segment if r_n > p_e
id_snow <- which (r_n <= p_e)
# p is used to generate snow load
p <- runif(length(id_snow))
# Final generation
g_s[id_snow] <- (B + 1.0 / A * (-log(-NS * log(1.0 - p_e + p_e * p))))^2*rooffactor(100,0.68,0.22)/
((Wind_velocity_50)^2*rooffactor(100,0.68,0.22))
g_s
mean(g_s)
##----------------------generate snow load-------------
N <- 100
# zip modified so that it can be applied to two lists
zip_modified <- function(a,b){
idx <- order(c(seq_along(a),seq_along(b)))
unlist(c(a,b)[idx])
}
# roof factor, convert ground snow load to roof snow load
# this is the parameter for vancouver
r_flat_mean <- 0.6
r_flat_COV <- 0.45
sigmalog <- function(mu,  sigma){
result <- sqrt(log(1 + exp(log(sigma*sigma) - 2 * log(mu))));
result
}
rooffactor <- function(N,r_flat_mean,r_flat_COV){
mu <- r_flat_mean
sigma <- r_flat_COV * r_flat_mean
input1 <- log(mu) - sigmalog(mu, sigma)*sigmalog(mu, sigma) / 2
input2 <- sigmalog(mu, sigma)
rlnorm(N,meanlog = input1, sdlog = input2)
}
# NS is the number of segment in each year
NS <- 10
# snow generating function
snow_generating<- function(N,NS){
# A = 4.3915  B = 0.1115 in vancouver
A <- 4.3915
B <- 0.1115
# p_0 is the probability of no-snow in one year (10 segment)
# p_e = 0.15 for vancouver
p0 <- exp(-exp(A*B))
# p_e is the probability of snow in a segment
p_e <- 1.00 - exp(-1.00/NS*(exp(A*B)))
# parameters for distribution g
Bstr <- A*B/(A*B+3.9019)
Astr <- A*B+3.9019
# generate N random number to decide whether it snows in one segment
r_n <- runif(N)
# generated snow laod
g_s <- rep(0,N)
# id that it snows in that segment. with r_n < p_e
# it does not snow in that winter segment if r_n > p_e
id_snow <- which (r_n <= p_e)
# p is used to generate snow load
p <- runif(length(id_snow))
# roof factor that convert the ground snow load to the roof snow load
r <- rooffactor(length(p),r_flat_mean = r_flat_mean, r_flat_COV = r_flat_COV)
# Final generation
g_s[id_snow] <- r*(Bstr + 1.0 / Astr * (-log(-NS * log(1.0 - p_e + p_e * p))))
g_s
}
##----------------------generate snow load-------------
N <- 100
# zip modified so that it can be applied to two lists
zip_modified <- function(a,b){
idx <- order(c(seq_along(a),seq_along(b)))
unlist(c(a,b)[idx])
}
# roof factor, convert ground snow load to roof snow load
# this is the parameter for vancouver
r_flat_mean <- 0.6
r_flat_COV <- 0.45
sigmalog <- function(mu,  sigma){
result <- sqrt(log(1 + exp(log(sigma*sigma) - 2 * log(mu))));
result
}
rooffactor <- function(N,r_flat_mean,r_flat_COV){
mu <- r_flat_mean
sigma <- r_flat_COV * r_flat_mean
input1 <- log(mu) - sigmalog(mu, sigma)*sigmalog(mu, sigma) / 2
input2 <- sigmalog(mu, sigma)
rlnorm(N,meanlog = input1, sdlog = input2)
}
# NS is the number of segment in each year
NS <- 10
# snow generating function
snow_generating<- function(N,NS){
# A = 4.3915  B = 0.1115 in vancouver
A <- 0.3222
B <- 17.0689
# p_0 is the probability of no-snow in one year (10 segment)
# p_e = 0.15 for vancouver
p0 <- exp(-exp(A*B))
# p_e is the probability of snow in a segment
p_e <- 1.00 - exp(-1.00/NS*(exp(A*B)))
# parameters for distribution g
Bstr <- A*B/(A*B+3.9019)
Astr <- A*B+3.9019
# generate N random number to decide whether it snows in one segment
r_n <- runif(N)
# generated snow laod
g_s <- rep(0,N)
# id that it snows in that segment. with r_n < p_e
# it does not snow in that winter segment if r_n > p_e
id_snow <- which (r_n <= p_e)
# p is used to generate snow load
p <- runif(length(id_snow))
# roof factor that convert the ground snow load to the roof snow load
r <- rooffactor(length(p),r_flat_mean = r_flat_mean, r_flat_COV = r_flat_COV)
# Final generation
g_s[id_snow] <- r*(Bstr + 1.0 / Astr * (-log(-NS * log(1.0 - p_e + p_e * p))))
g_s
}
T_snow_between <- rexp(N, 1/(7/12))
# generate the duration time (two weeks per segment)
# generate the snow load of each segment in the winter
require(tidyverse)
T_snow_duration <- vector(mode = 'list',length = N) %>%
lapply(function(x) x = rexp(NS, 1/0.03835))
load_snow <- vector(mode = 'list',length = N) %>%
lapply(function(x) x = snow_generating(NS))
# generate the final snow laod
D_snow <- stepfun(cumsum(zip_modified(T_snow_between, T_snow_duration)),
zip_modified(rep(0, N+1), load_snow))
snow_generating(N = 100, NS =10)
mean(snow_generating(N = 100, NS =10))
if (require("remotes", quietly = TRUE) == FALSE) {
install.packages("remotes")
require("remotes")
}
remotes::install_github("vdorie/aciccomp/2016")
load("~/Downloads/theta/GP-llikObs.rda")
colMeans(theta)
load("~/Downloads/github/damage-models/fit-gp/1.rda")
colMeans(theta[1:100,])
View(theta)
load("~/Downloads/theta/GP-llikObs.rda")
load("~/Downloads/github/damage-models/fit-gp/1.rda")
View(theta)
colMeans(theta)
colMeans(theta[1:1200,])
load("~/Downloads/theta/GP-llikObs.rda")
colMeans(theta)
colMeans(theta)[2]
colMeans(theta)[,3]
colMeans(theta)[3]
install.packages("qte")
library(qte)
data(lalonde)
force(lalonde.psid.panel)
install.packages("causalsens")
install.packages("coda")
install.packages("bayesGARCH")
install.packages("coda")
require(bayesGARCH)
## LOAD DATA
data(dem2gbp)
y <- dem2gbp[1:750]
## RUN THE SAMPLER (2 chains)
MCMC <- bayesGARCH(y, control = list(n.chain = 2, l.chain = 200))
## MCMC ANALYSIS (using coda)
plot(MCMC)
par(mar = rep(2, 4))
plot(MCMC)
## LOAD DATA
data(dem2gbp)
y <- dem2gbp[1:750]
## RUN THE SAMPLER (2 chains)
par(mar = rep(2, 4))
MCMC <- bayesGARCH(y, control = list(n.chain = 2, l.chain = 500))
## MCMC ANALYSIS (using coda)
plot(MCMC)
## LOAD DATA
data(dem2gbp)
y <- dem2gbp[1:750]
## RUN THE SAMPLER (2 chains)
par(mar = rep(2, 4))
MCMC <- bayesGARCH(y, control = list(n.chain = 4, l.chain = 500))
## MCMC ANALYSIS (using coda)
plot(MCMC)
windows() ## create window to plot your file
plot(MCMC)
dev.off()
windows() ## create window to plot your file
Sys.setenv(LANGUAGE = "en")
windows() ## create window to plot your file
window() ## create window to plot your file
dev.new() ## create window to plot your file
plot(MCMC)
dev.off()
## LOAD DATA
data(dem2gbp)
y <- dem2gbp[1:750]
## RUN THE SAMPLER (2 chains)
MCMC <- bayesGARCH(y, control = list(n.chain = 4, l.chain = 500),
lambda = 100,delta = 500)
## MCMC ANALYSIS (using coda)
par(mar = rep(2, 4))
plot(MCMC)
View(MCMC)
View(MCMC)
MCMC[1]
draw <- MCMC$chain1
ts.plot(draws[range,1],xlab="iterations",main=expression(alpha[0]),ylab="")
ts.plot(draw[range,1],xlab="iterations",main=expression(alpha[0]),ylab="")
ts.plot(draw[,1],xlab="iterations",main=expression(alpha[0]),ylab="")
plot(dem2gbp[1:750])
plot(dem2gbp[1:750],type = 'l')
getwd()
source('R/fit_US.R
')
source("R/fit_US.R")
fit1s
rooffactor(100,0.68,0.22)
?quantile
a <- rooffactor(100000,0.68,0.22)
quantile(a,probs = 0.95)
quantile(a,probs = 0.90)
quantile(a,probs = 49/50)
mean(a)
G_cov <- 0.150
G_mean <- (1+3.050*G_cov)/(1+2.592*G_cov)
r_flat_mean = G_mean
r_flat_COV = G_cov
mu <- r_flat_mean
sigma <- r_flat_COV * r_flat_mean
input1 <- log(mu) - sigmalog(mu, sigma)*sigmalog(mu, sigma) / 2
input2 <- sigmalog(mu, sigma)
exp(input1 + input2^2/2)
r_flat_mean = 0.68
r_flat_COV = 0.22
mu <- r_flat_mean
sigma <- r_flat_COV * r_flat_mean
input1 <- log(mu) - sigmalog(mu, sigma)*sigmalog(mu, sigma) / 2
input2 <- sigmalog(mu, sigma)
exp(input1 + input2^2/2)
load("~/Downloads/GP-llikObs.rda")
remove(llikObs)
setwd("~/Downloads")
save.image("~/Downloads/GP-theta.RData")
save.image("~/Downloads/GP-theta.RData",version = 2)
save.image("~/Downloads/GP-theta.RData",version = 2)
save(theta,version = 2,file = 'GP-theta.rda')
require(bayesGARCH)
data("dem2gbp")
?dem2gbp
plot(dem2gbp,type = 'l')
plot(dem2gbp,type = 'l',main = 'Foreign exchange rate')
plot(dem2gbp[1:1200],type = 'l',main = 'Foreign exchange rate')
plot(dem2gbp,type = 'l',main = 'Foreign exchange rate')
setwd("~/Downloads/github/Multimodel-Reliability")
tau0 <- 6468 # Median short-term strength
##  Setup
source("R/loadData.R")
numData <- length(obs)
source("R/bunches.R")
EOS <- list()
for (i in 1:numData) {
EOS[[i]] <- EnvStats::evNormOrdStats(length(obs[[i]]))
}
wNLS <- function(logT, R, k, tauc, Tp0, Tp1, case, wgt, a, b, w) {
ewR <- exp(w * R)
alpha_T1 <- ifelse(case ==2, (Tp1 - Tp0) * exp ( -a + b * tauc / ewR) + exp(-a) * ewR / (b*k) * ( exp( b * k * Tp0 / ewR) - 1), 0)
res <- ifelse( case == 0, logT - log(ewR / (b * k) * log ( exp(a) * b * k / ewR + 1)),
ifelse ( case == 1, (logT - log(tauc/k - ewR / (b*k) + exp( -b * tauc / ewR) * ( ewR / (b*k) + exp(a)))) * wgt,
logT - log( ewR / (b * k) * log ( exp(a) * b * k / ewR * (1 - alpha_T1)+ 1)   ) )   #case = 2
)
if(sum(alpha_T1>1)>0) message(paste("alpha_T1 exceeded in sample# ", which(alpha_T1>1)))  ## bad cases
ifelse( is.nan(res), (alpha_T1-1)*1000, res)
}
alldat <- data.frame(T = (unlist(obs[1:numData])),
R = unlist(EOS[1:numData]),
k = unlist(sapply(1:numData, function(i) rep(k_app[i], length(obs[[i]])))),
i = unlist(sapply(1:numData, function(x) rep(x, length(obs[[x]])))),
tauc = unlist(sapply(1:numData, function(x) rep(tauc[x], length(obs[[x]])))),
Tp0 = unlist(sapply(1:numData, function(x) rep(T0[x], length(obs[[x]])))),
Tp1 = unlist(sapply(1:numData, function(x) rep(T1[x], length(obs[[x]])))),
wgt = 1
)
alldat$logT <- with(alldat, {ifelse( T > T1[i] & T1[i] > 0, log(T - T1[i]), log(T) )})
alldat$case <- with(alldat, {ifelse( T > T1[i] & T1[i] > 0, 2, ifelse(T > T0[i] & T0[i] > 0, 1, 0))})
fit1 <- nls( ~ wNLS(logT, R, k, tauc, Tp0, Tp1, case, wgt, a, b, w), data = alldat, start=list(a=20.665529538, b = 0.000183925, w = 0.803524486), trace=T, control=list(maxiter=100,printEval = T))
fit1s <- summary(fit1)
itcursigma <- 1e9
alldat$wgt <- ifelse(alldat$case==1, 1/fit1s$coefficients[2]/alldat$tauc , 1)
### Iterate until convergence
while( abs(itcursigma - fit1s$sigma) > 1e-8) {
itcursigma <- fit1s$sigma
fit1 <- nls( ~ wNLS(logT, R, k, tauc, Tp0, Tp1, case, wgt, a, b, w), data = alldat, start=list(a = fit1s$coefficients[1], b = fit1s$coefficients[2], w = fit1s$coefficients[3]), trace=T, control=list(maxiter=1000,printEval = T))
fit1s <- summary(fit1)
alldat$wgt <- ifelse(alldat$case==1, 1/fit1s$coefficients[2]/alldat$tauc , 1)
}
ger_a <- fit1s$coefficients[1]
ger_b <- fit1s$coefficients[2]  # this is actually B' (we use it throughout in USADM as well)
ger_w <- fit1s$coefficients[3]
## Write out parameter estimate
write.table( cbind(ger_a,ger_b,ger_w), file="US_theta.csv", row.names = F, sep=',', col.names = F)
fit1s$coefficients
load("~/Downloads/US-llikObs.rda")
colMeans(theta)
load("~/Downloads/theta/Can-llikObs.rda")
theta[1,]
setwd()
getwd()
write.table( thetam file="Can_theta.csv", row.names = F, sep=',', col.names = F)
write.table( theta, file="Can_theta.csv", row.names = F, sep=',', col.names = F)
library(deSolve)
# Number of stochastic load profiles to generate
# nrep <- 100000
nrep <- 20
# Load parameters
gamma = 0.25; R_0 = 3000; alpha_d = 1.25; alpha_l = 1.5
## tau(t):  this will be stochastically generated for live load applications
zip <- function(a, b){
idx <- order(c(seq_along(a), seq_along(b)))
unlist(c(a,b))[idx]
}
D_d <- list()
T_s <- list()
load_s <- list()
D_s <- list()
T_e <- list()
T_p <- list()
load_p <- list()
D_e <- list()
N <- 100
for (i in 1:nrep) {
D_d[[i]] <- rnorm(1, 1.05, 0.1)
T_s[[i]] <- rexp(N, 0.1)
load_s[[i]] <- rgamma(N+1, shape = 3.122, scale = 0.0481)
D_s[[i]] <- stepfun(cumsum(T_s[[i]]), load_s[[i]])
T_e[[i]] <- rexp(N, 1)
T_p[[i]] <- rexp(N, 1/0.03835)
load_p[[i]] <- rgamma(N, shape = 0.826, scale = 0.1023)
D_e[[i]] <- stepfun(cumsum(zip(T_e[[i]], T_p[[i]])), zip(rep(0, N+1), load_p[[i]]))
}
length(zip(rep(0, N+1), load_p[[1]]))