Practice Histograms

#Random Normal Histogram
random <- rnorm(100, mean = 0, sd = 1)
hist(random, freq = FALSE, main = "Histogram of Normal Fitted Curve")
x<- seq(min(random), max(random), length = 100)
sdensity <- dnorm(x, mean = 0, sd = 1, log = FALSE)
lines(x, sdensity, col = "red")
mle <- dnorm(x, mean = mean(random), sd = sd(random), log = FALSE)
lines(x, mle, col = "blue")

#Random Exponential Histogram
erandom <- rexp(100, rate = 1)
hist(erandom, freq = FALSE, main = "Histogram of Exponential Fitted Curve")
x<- seq(min(erandom), max(erandom), length = 100)
edensity <- dexp(x, rate = 1, log = FALSE)
lines(x, edensity, col = "red")
mle <- dexp(x, rate = 1/mean(erandom), log = FALSE)
lines(x, mle, col = "blue")

Plots and Histograms

#PrX <= x = F(x) = 1 - [1 + shape * (x - threshold)/scale]^(-1/shape)
#F(x) = 1 - exp(-(x - threshold)/scale)

#Differing Histograms with Exponential and Pareto MLE
STATION <- "SAN LUIS OBISPO POLY CA US"    #enter/change station name here
PRCP <- precip$SAN.LUIS.OBISPO.POLY.CA.US  #enter/change station name here
PRCP <- subset(PRCP, PRCP != "NA")         #removing NAs
PRCP <- as.numeric(PRCP)                   #converting from char strings to numbers
SPRCP <- PRCP[PRCP >= 0.01]                
xgrid<- seq(min(SPRCP), max(SPRCP), length = length(SPRCP))
xgrid <- xgrid[xgrid >= 0.01]
fit <- fevd(PRCP, threshold = 0.1, type="GP")
title <- paste("Histogram of ", STATION, " Precipitation with Fitted Curves")
hist(SPRCP, freq = FALSE, main = title, xlab = "Precipitation (in)", ylab = "Density")
dens <- devd(xgrid, scale = fit$results$par[1], shape = fit$results$par[2], threshold = fit$threshold, log = FALSE, type = "GP")
mle <- dexp(xgrid, rate = 1/mean(SPRCP), log = FALSE)
lines(xgrid, mle, col = "blue")            #exponential fitted curve
lines(xgrid, dens, col = "red")            #Pareto fitted curve
legend("topright", c("Exponential", "Generalized Pareto"), lty = c(1,1), col = c("blue", "red"))

#Differing Return Levels for 100 Year Return Period per station
p <- precip$SAN.LUIS.OBISPO.POLY.CA.US  #enter/change station name here
p <- subset(p, p != "NA")
p <- as.numeric(p)
fit <- fevd(p, threshold = 1, type="GP")
fit
## 
## fevd(x = p, threshold = 1, type = "GP")
## 
## [1] "Estimation Method used: MLE"
## 
## 
##  Negative Log-Likelihood Value:  469.8975 
## 
## 
##  Estimated parameters:
##      scale      shape 
## 0.72684315 0.01313252 
## 
##  Standard Error Estimates:
##      scale      shape 
## 0.03991997 0.03923778 
## 
##  Estimated parameter covariance matrix.
##              scale        shape
## scale  0.001593604 -0.001104778
## shape -0.001104778  0.001539604
## 
##  AIC = 943.7951 
## 
##  BIC = 952.8304
plot(fit)

fit$results$par
##      scale      shape 
## 0.72684315 0.01313252
pthres <- sum(p > 1)/length(p) 
qevd(1 / (pthres*(100*365)), scale = fit$results$par[1],      #Return Level for 100 Years
     shape = fit$results$par[2], threshold = fit$threshold,
       type = "GP", lower.tail = FALSE)
##    scale 
## 6.012011

Return Level, Scale, and Shape Plots on Map

#Locations
new_locations <- subset(locations, lon != "unknown") #one station has no location, removed from data set
station_ids2 <- row.names(new_locations)
#maptype options are “terrain”, “satellite”, “roadmap”, and “hybrid”
map2 <- get_map(location = c(-120.5, 35.4), zoom = 9, maptype = "terrain", scale = 2)
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=35.4,-120.5&zoom=9&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
#scalecolors
##Return Levels with threshold = 1
ggmap(map2) + geom_point(alpha = 1, aes(Longitude, Latitude, color = Return.Level, size = Size), data = thres1) + labs(title =  "San Luis Obispo County Station Return Levels w/ Threshold 1", x = "Longitude", y = "Latitude") + scale_colour_gradient(low = "blue", high = "red")

##Scales with threshold = 1
ggmap(map2) + geom_point(alpha = 1, aes(Longitude, Latitude, color = Scale, size = Size), data = thres1) + labs(title =  "San Luis Obispo County Station Return Levels w/ Threshold 1", x = "Longitude", y = "Latitude") + scale_colour_gradient(low = "blue", high = "red")

##Shapes with threshold = 1
ggmap(map2) + geom_point(alpha = 1, aes(Longitude, Latitude, color = Shape, size = Size), data = thres1) + labs(title =  "San Luis Obispo County Station Return Levels w/ Threshold 1", x = "Longitude", y = "Latitude") + scale_colour_gradient(low = "blue", high = "red")

#Locations
new_locations <- subset(locations, lon != "unknown") #one station has no location, removed from data set
station_ids2 <- row.names(new_locations)
#maptype options are “terrain”, “satellite”, “roadmap”, and “hybrid”
map2 <- get_map(location = c(-120.5, 35.4), zoom = 9, maptype = "terrain", scale = 2)
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=35.4,-120.5&zoom=9&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
#rainbowcolors
##Return Levels with threhold = 1
ggmap(map2) + geom_point(alpha = 1, aes(Longitude, Latitude, color = factor(Return.Level), size = Size), data = thres1) + labs(title =  "San Luis Obispo County Station Return Levels w/ Threshold 1", x = "Longitude", y = "Latitude")

##Scales with threhold = 1
ggmap(map2) + geom_point(alpha = 1, aes(Longitude, Latitude, color = factor(Scale), size = Size), data = thres1) + labs(title =  "San Luis Obispo County Station Return Levels w/ Threshold 1", x = "Longitude", y = "Latitude")

##Shapes with threhold = 1
ggmap(map2) + geom_point(alpha = 1, aes(Longitude, Latitude, color = factor(Shape), size = Size), data = thres1) + labs(title =  "San Luis Obispo County Station Return Levels w/ Threshold 1", x = "Longitude", y = "Latitude")