Weibull distribution and wind data.
winddata <-read.csv("plotweibull_excercise.csv",
header=TRUE, sep=",")
str(density(winddata$X10m))
## List of 7
## $ x : num [1:512] -1.92 -1.84 -1.75 -1.67 -1.59 ...
## $ y : num [1:512] 1.50e-05 2.21e-05 3.20e-05 4.56e-05 6.38e-05 ...
## $ bw : num 0.639
## $ n : int 32512
## $ call : language density.default(x = winddata$X10m)
## $ data.name: chr "winddata$X10m"
## $ has.na : logi FALSE
## - attr(*, "class")= chr "density"
# extract the density information
my.x <- density(winddata$X10m, from = 0, to = max(winddata$X10m))$x
my.y <- density(winddata$X10m, from = 0, to = max(winddata$X10m))$y
# assemble a data frame
df <- data.frame(x = my.x, y = my.y)
plot(y ~ x, data = df)
# plot with additional information,
# this plot is shown at the end
plot(density(winddata$X10m, from = 0, to = max(winddata$X10m)))
# adding the weibull distribution, with information on shape and scape from a previous analysis
curve(dweibull(x, scale = 15.47, shape = 2.37), from = 0, to = 45, add = TRUE, col = "red")
# manual curve fit via nls
# shape = a, scale = b
# well, weibull works only for x > 0!
my.nls <- nls(y ~ (a/b) * (x/b)^(a-1) * exp(- (x/b)^a),
data = df[df$x > 0, ],
start = c(a = 2.37, b = 15.47))
summary(my.nls)
##
## Formula: y ~ (a/b) * (x/b)^(a - 1) * exp(-(x/b)^a)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 2.6000 0.0142 183 <2e-16 ***
## b 15.0369 0.0408 369 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.00387 on 509 degrees of freedom
##
## Number of iterations to convergence: 8
## Achieved convergence tolerance: 2.56e-06
xValues <- seq(from = 0, to = 40, length.out = 100) # create a sequence from 0 to the desired max x-value with 100 elements
# predict values from the nls.fit.
# Predict expects a dataframe "in which to look for variables with which to predict".
# Name has to match the nls model for the independent variable / predictor.
my.predicted <- predict(my.nls, data.frame(x = xValues))
# use lines to connect the predicted values.
lines(xValues,
my.predicted,
col = "blue")
# grab the coefficients
my.coef <- coef(my.nls)
# a little function that gives the y-value back for a given x and given coefficients.
my.weibull.predict <- function(x, a, b) {
y <- (a/b) * (x/b)^(a-1) * exp(- (x/b)^a)
return(y)
}
# add some points to the plot
my10 <- my.weibull.predict(10, my.coef[1], my.coef[2])
my29 <- my.weibull.predict(29, my.coef[1], my.coef[2])
points(10, my10, col = "green")
points(29, my29, col = "red")
