Moyenne d'une fenêtre coulissante en R

19

J'ai un vecteur de valeurs que je voudrais signaler la moyenne dans les fenêtres le long d'une petite diapositive.

Par exemple, pour un vecteur des valeurs suivantes:

4, 5, 7, 3, 9, 8

Une taille de fenêtre de 3 et une diapositive de 2 feraient ce qui suit:

(4+5+7)/3 = 5.33
(7+3+9)/3 = 6.33
(9+8)/3 = 5.67

Et retournez un vecteur de ces valeurs:

5.33, 6.33, 5.67

Y a-t-il une fonction simple qui fera cela pour moi? S'il renvoyait également les indices de la fenêtre, ce serait un bonus supplémentaire. Dans cet exemple, ce serait 1,3,5

T-Burns
la source
4
Tu as vu ça ?
JM n'est pas statisticien le
Pouvez-vous donner un aperçu de cette idée de "diapositive"?
Shane
@JM - Je ne l'avais pas fait! Je vous remercie! Je vais voir comment ça marche.
T-Burns du
@Shane - Oui! Je suis désolé, ce n'était pas clair. La diapositive est le nombre de positions / indices que vous déplacez pour commencer à calculer la prochaine fenêtre de moyennes. Ainsi, plutôt que la fenêtre suivante commençant après la fin de la dernière, il y a un certain chevauchement lorsque la diapositive est plus petite que la taille de votre fenêtre. L'idée est de lisser un peu les points de données.
T-Burns
Merci, j'avais la même question. Maintenant, j'ai trouvé utile la fonction "rollapply".
angelous

Réponses:

24

La fonction rollapplydu package zoo vous rapproche:

> require(zoo)
> TS <- zoo(c(4, 5, 7, 3, 9, 8))
> rollapply(TS, width = 3, by = 2, FUN = mean, align = "left")
       1        3 
5.333333 6.333333

Il ne calculera tout simplement pas la dernière valeur pour vous car il ne contient pas 3 observations. Peut-être que ce sera suffisant pour votre vrai problème? Notez également que l'objet renvoyé a les indices que vous souhaitez en tant que du namesvecteur renvoyé.

Votre exemple fait l'hypothèse qu'il y a un 0 non observé dans la dernière fenêtre. Il peut être plus utile ou réaliste de remplir avec un NApour représenter les informations manquantes et de dire meande gérer les valeurs manquantes. Dans ce cas, nous aurons (8 + 9) / 2 comme valeur fenêtrée finale.

> TS <- zoo(c(4, 5, 7, 3, 9, 8, NA))
> rollapply(TS, width = 3, by = 2, FUN = mean, na.rm = TRUE, align = "left")
       1        3        5 
5.333333 6.333333 8.500000
Réintégrer Monica - G. Simpson
la source
BTW, j'ai écrit une fois sur l'utilisation de cette fonction pour implémenter la notion de "loess quantile": r-statistics.com/2010/04/…
Tal Galili
Vous pouvez ajouter un 0 à la fin de x ( x<-c(x,0)) pour obtenir le dernier élément de réponse.
1
@mbq; cela fait une forte hypothèse que l'observation est 0. J'avais réfléchi à ce point et T-Burns fait la même hypothèse (un 0 non observé). Je préférerais peut-être jouer avec NA et passer l' na.rm = TRUEargument à mean. La réponse ne sera pas la même que celle demandée par le PO, mais elle semble plus utile. Je vais modifier ma réponse pour l'inclure.
Rétablir Monica - G. Simpson
@ucfagls Pourtant, cela est facile à changer et comme vous l'avez dit, cette hypothèse a été émise par le PO. En revanche, je serais encore plus restrictif et supprimerais la dernière moyenne.
Merci! Surtout pour avoir noté la dernière valeur comme une hypothèse nulle, je n'avais pas considéré cela. Je me soucie vraiment de cette dernière fenêtre !!
T-Burns du
12

Rollapply fonctionne très bien avec un petit ensemble de données. Cependant, si vous travaillez avec plusieurs millions de lignes (génomique), c'est assez lent.

La fonction suivante est super rapide.

data <- c(runif(100000, min=0, max=.1),runif(100000, min=.05, max=.1),runif(10000, min=.05, max=1), runif(100000, min=0, max=.2))

slideFunct <- function(data, window, step){
  total <- length(data)
  spots <- seq(from=1, to=(total-window), by=step)
  result <- vector(length = length(spots))
  for(i in 1:length(spots)){
    result[i] <- mean(data[spots[i]:(spots[i]+window)])
  }
  return(result)
}

http://coleoguy.blogspot.com/2014/04/sliding-window-analysis.html

r_evolutionist
la source
Très utile. Mais sachez que cette fenêtre = 3 renverra la moyenne de 4 (!) Valeurs, sauf si vous ajoutez un -1(à la plage) et un +1(à la boucle).
BurninLeo
5

Cette simple ligne de code fait la chose:

((c(x,0,0) + c(0,x,0) + c(0,0,x))/3)[3:(length(x)-1)]

si xest le vecteur en question.

user1414
la source
Cela ne renvoie pas ce que le demandeur voulait, mais 5,33 5,00 6,33. Cependant, cela semble assez intéressant. Pouvez-vous expliquer votre idée, car je ne comprends pas.
Henrik
1
@Henric J'utilise cette astuce fréquemment, mais le code de user1414 renvoie ce rouleau avec la diapositive 1, pas 2, comme prévu par OP. Découvrez (c(0,0,x)+c(0,x,0)+c(x,0,0))/3ce que je veux dire (et comment ça marche). La formule appropriée serait: (c(0,0,x)+c(0,x,0)+c(x,0,0))[1:(length(x)-3)*2+1]/3(nous devons couper le remplissage 0 au début et sélectionner ensuite les éléments pairs.
4
library(zoo)
x=c(4, 5, 7, 3, 9, 8)
rollmean(x,3)

ou

library(TTR)
x=c(4, 5, 7, 3, 9, 8)
SMA(x,3)
RockScience
la source
Est-ce que cela fonctionne pour les matrices 2D? Comme quoi? Si la taille de la fenêtre est de 3 * 3 à titre d'exemple
Mona Jalal
ce n'est qu'une seule direction
RockScience
3

réponse de shabbychef dans R:

slideMean<-function(x,windowsize=3,slide=2){
 idx1<-seq(1,length(x),by=slide);
 idx1+windowsize->idx2;
 idx2[idx2>(length(x)+1)]<-length(x)+1;
 c(0,cumsum(x))->cx;
 return((cx[idx2]-cx[idx1])/windowsize);
}

EDIT: Les indices que vous recherchez sont juste idx1... cette fonction peut être facilement modifiée pour les renvoyer également, mais il est presque aussi rapide de les recréer avec un autre appel à seq(1,length(x),by=slide).

Communauté
la source
merci d'avoir traduit. Je pensais que ce serait un exercice facile, et j'en ai appris du R
shabbychef
Ma réponse mise à jour est l'utilisation fromo::running_meande la version de pointe de mon package fromo .
shabbychef
3

Je peux le faire facilement dans Matlab et canard pendant que vous me downvote:

%given vector x, windowsize, slide 
idx1 = 1:slide:numel(x);
idx2 = min(numel(x) + 1,idx1 + windowsize);  %sic on +1 here and no -1;
cx = [0;cumsum(x(:))];  %pad out a zero, perform a cumulative sum;
rv = (cx(idx2) - cx(idx1)) / windowsize; %tada! the answer!

comme effet secondaire, idx1est l'indice de l'élément dans la somme. Je suis sûr que cela peut être facilement traduit en R. L'idiome first:skip:lastdans Matlab donne le tableau en premier, premier + saut, premier + 2 saut, ..., premier + n saut, où le dernier élément du tableau n'est pas supérieur àlast .

edit : j'avais omis la partie moyenne (diviser par windowsize).

shabbychef
la source
+1 Pas tada, rv /
1
Cette zone de commentaire marg ... est trop étroite pour ce code, j'ai donc posté une nouvelle réponse.
1
Merci, mais MATLAB n'est pas gratuit !!
T-Burns du
@ T-Burns: l'octave est cependant libre; R est également assez proche de Matlab pour que ce code puisse être facilement traduit. En fait, @mbq a fait ça ..
shabbychef
1

Cela vous donnera les moyennes de la fenêtre et l'index de la première valeur de la fenêtre:

#The data
x <- c(4, 5, 7, 3, 9, 8)

#Set window size and slide
win.size <- 3
slide <- 2

#Set up the table of results
results <- data.frame(index = numeric(), win.mean = numeric())

#i indexes the first value of the window (the sill?)
i <- 1
#j indexes the row of the results to be added next
j <- 1
while(i < length(x)) {
    #This mean preserves the denominator of 3
    win.mean <- sum(x[i:(i+2)], na.rm = TRUE)/win.size
    #Insert the results
    results[j, ] <- c(i, win.mean)
    #Increment the indices for the next pass
    i <- i + slide
    j <- j + 1
    }

Diverses mises en garde s'appliquent: je n'ai pas testé cela par rapport à vos données d'échantillon; Je crois que l'ajout à des trames de données comme celle-ci peut devenir très lent si vous avez beaucoup de valeurs (car cela copiera le data.frame à chaque fois); etc. Mais cela produit ce que vous avez demandé.

Matt Parker
la source
S'il vous plaît, ne réduisez pas votre note sans fournir de commentaire. Comment suis-je censé savoir ce qui ne va pas?
Matt Parker
Ce n'était pas moi, mais c'est lent (mais pas beaucoup plus lent que rollapply).
2
ce n'était pas moi non plus, mais comme vous l'avez mentionné, la pré-allocation de l'objet résultat aidera à résoudre le problème de vitesse. Une astuce, si vous ne connaissez pas, ou si c'est fastidieux / difficile à déterminer, la taille de l'objet de résultat dont vous avez besoin. Allouez quelque chose de raisonnable, peut-être en pré-remplissant NA. Remplissez ensuite votre boucle, mais ajoutez une vérification que si vous approchez de la limite de l'objet préalloué, allouez un autre gros morceau et continuez à remplir.
Rétablir Monica - G. Simpson
1
@mbq; La rapidité des résultats, bien qu'importante, n'est pas la seule considération. Au lieu d'avoir à réinventer le temps et à gérer tous les index, etc. dans les solutions personnalisées, celui linéaire rollapplyest beaucoup plus facile à comprendre et à comprendre l'intention de. En outre, il rollapplyest probable qu'il y ait eu beaucoup plus de globes oculaires vérifiant son code que quelque chose que je pourrais cuisiner un après-midi. Chevaux de course.
Rétablir Monica - G. Simpson,
1
Changer [i:(i+2)]pour [i:(i+win.size-1)]rendre le code plus général, je pense.
Jota