Добавление сезонных колебаний во временные ряды скорости ветра

Question

В продолжение блога R, который интересен и весьма полезен для моделирования временных рядов неизвестной области с использованием ее параметров Вейбулла.

Хотя этот метод дает достаточно хорошую оценку временных рядов в целом, он сильно страдает, когда мы ищем сезонные изменения. Чтобы учесть сезонные изменения, я хочу использовать сезонные максимальные скорости ветра и провести синтез временных рядов так, чтобы годовое распределение оставалось постоянным, т.е. параметры формы и масштаба (годовые значения).

Я хочу использовать сезонные максимальные скорости ветра в приведенном ниже коде, используя 12 различных максимальных скоростей ветра, по одной на каждый месяц. Это позволит увеличить скорость ветра в определенный месяц и снизить в другие, а также выровнять результирующий временной ряд.

Код выглядит следующим образом:

MeanSpeed<-7.29 ## Mean Yearly Wind Speed at the site.

Shape=2; ## Input Shape parameter (yearly).
Scale=8 ##Calculated Scale Parameter ( yearly).

MaxSpeed<-17 (##yearly)
## $$$ 12 values of these wind speed one for each month to be used. The resultant time series should satisfy shape and scale parameters $$ ###
nStates<-16 

nRows<-nStates;
nColumns<-nStates;


LCateg<-MaxSpeed/nStates; 


WindSpeed=seq(LCateg/2,MaxSpeed-LCateg/2,by=LCateg) ## Fine the velocity vector-centered on the average value of each category.

##Determine Weibull Probability Distribution.
wpdWind<-dweibull(WindSpeed,shape=Shape, scale=Scale); # Freqency distribution.

plot(wpdWind,type = "b", ylab= "frequency", xlab = "Wind Speed")  ##Plot weibull probability distribution.

norm_wpdWind<-wpdWind/sum(wpdWind); ## Convert weibull/Gaussian distribution to normal distribution.

## Correlation between states (Matrix G)
g<-function(x){2^(-abs(x))} ## decreasing correlation function between states.
G<-matrix(nrow=nRows,ncol=nColumns)
G <- row(G)-col(G)
G <- g(G)

##--------------------------------------------------------


## iterative process to calculate the matrix P (initial probability)
P0<-diag(norm_wpdWind);   ## Initial value of the MATRIX P.
P1<-norm_wpdWind;  ## Initial value of the VECTOR p.


## This iterative calculation must be done until a certain error is exceeded
## Now, as something tentative, I set the number of iterations

steps=1000;  
P=P0; 
p=P1; 

for (i in 1:steps){
    r<-P%*%G%*%p;
    r<-as.vector(r/sum(r)); ## The above result is in matrix form. I change it to vector
    p=p+0.5*(P1-r)
    P=diag(p)}

   ## $$ ----Markov Transition Matrix --- $$ ##

N=diag(1/as.vector(p%*%G));## normalization matrix

MTM=N%*%G%*%P ## Markov Transition Matrix

MTMcum<-t(apply(MTM,1,cumsum));## From the MTM generated the accumulated

##-------------------------------------------
## Calculating the series from the MTMcum

##Insert number of data sets. 
LSerie<-52560; Wind Speed every 10 minutes for a year. 

RandNum1<-runif(LSerie);## Random number to choose between states
State<-InitialState<-1;## assumes that the initial state is 1 (this must be changed when concatenating days)
StatesSeries=InitialState;

## Initallise----

## The next state is selected to the one in which the random number exceeds the accumulated probability value
##The next iterative procedure chooses the next state whose random number is greater than the cumulated probability defined by the MTM
for (i in 2:LSerie) {
  ## i has to start on 2 !!
    State=min(which(RandNum1[i]<=MTMcum[State,]));

    ## if (is.infinite (State)) {State = 1}; ## when the above condition is not met max -Inf
    StatesSeries=c(StatesSeries,State)}

RandNum2<-runif(LSerie); ## Random number to choose between speeds within a state

SpeedSeries=WindSpeed[StatesSeries]-0.5+RandNum2*LCateg;
##where the 0.5 correction is needed since the the WindSpeed vector is centered around the mean value of each category.


print(fitdistr(SpeedSeries, 'weibull')) ##MLE fitting of SpeedSeries

Кто-нибудь может подсказать, где и какие изменения мне нужно внести в код?

Answer 1

Я не знаю много о генерации временных рядов скорости ветра, но, возможно, эти рекомендации помогут вам улучшить читаемость / повторное использование кода:

# 1 Вы, вероятно, хотите иметь функцию, которая будет генерировать серию времени скорости ветра, учитывая количество наблюдений и сезонную максимальную скорость ветра. Итак, сначала попробуйте определить свой код внутри блока, как этот:

wind_time_serie <- function(nobs, max_speed){
    #some code here
}

# 2 При этом, если кажется, что некоторые части вашего кода полезны для генерации временных рядов скорости ветра, но не относятся к временным рядам скорости ветра, попробуйте поместить их в функции (например, ту часть, которую вы вычисляете norm_wpdWindчасть, которую вы вычисляете MTMcum,...).

# 3 Затем часть вашего кода в начале, когда ваша определяемая глобальная переменная должна исчезнуть и стать аргументами по умолчанию в функциях.

# 4 Избегайте использования комментариев конца строки, когда ваша строка уже длинная, и удалите конечные точки с запятой.

#This
  State<-InitialState<-1;## assumes that the initial state is 1 (this must be changed when concatenating days)

#Would become this:
  #Assumes that the initial state is 1 (this must be changed when concatenating days)
  State<-InitialState<-1

Тогда ваш код должен быть более пригодным для повторного использования / чтения другими людьми. Ниже приведен пример тех рекомендаций, которые применяются к части rnorm:

  norm_distrib<-function(maxSpeed, states = 16, shape = 2, scale = 8){

    #Fine the velocity vector-centered on the average value of each category.
    LCateg<-maxSpeed/states
    WindSpeed=seq(LCateg/2,maxSpeed-LCateg/2,by=LCateg) 

    #Determine Weibull Probability Distribution.
    wpdWind<-dweibull(WindSpeed,shape=shape, scale=scale)

    #Convert weibull/Gaussian distribution to normal distribution.
    return(wpdWind/sum(wpdWind))

  }

  #Plot normal distribution with the max speed you want (e.g. 17)
  plot(norm_distrib(17),type = "b", ylab= "frequency", xlab = "Wind Speed")