CODEWITHSUNDEEP
r
balin
balin

Reputation: 1686

R: Finding the begin of a (exponential?) decay?

How to find the index indicated by the red vlin in the following example:

# Get the data as "tmpData"
source("http://pastie.org/pastes/9350691/download")

# Plot 
plot(tmpData,type="l")
abline(v=49,col="red")

First Plot

The following approach is promising, but how to find the peak maximum?

library(RcppRoll)
n <- 10
smoothedTmpData <- roll_mean(tmpData,n)
plot(-diff(smoothedTmpData),type="l")
abline(v=49,col="red")

Second Plot

Upvotes: 2

Views: 113

Answers (2)

gwieshammer
gwieshammer

Reputation: 21

which.max(-diff(smoothedTmpData)) gives you the index of the maximum.

http://www.inside-r.org/r-doc/base/which.max

I'm unsure if this is your actual question...

Upvotes: 2

Richie Cotton
Richie Cotton

Reputation: 121127

Where there is a single peak in the gradient, as in your example dataset, then gwieshammer is correct: you can just use which.max to find it.

For the case where there are multiple possible peaks, you need a more sophisticated approach. R has lots of peak finding functions (of varying quality). One that works for this data is wavCWTPeaks in wmtsa.

library(RcppRoll)
library(wmtsa)

source("http://pastie.org/pastes/9350691/download")

n <- 10
smoothedTmpData <- roll_mean(tmpData, n)

gradient <- -diff(smoothedTmpData)

cwt <- wavCWT(gradient)
tree <- wavCWTTree(cwt)
(peaks <- wavCWTPeaks(tree))
## $x
## [1]  4 52
## 
## $y
## [1]  302.6718 5844.3172
## 
## attr(,"peaks")
## branch itime iscale time scale  extrema iendtime
## 1      1     5      2    5     2 16620.58        4
## 2      2    57     26   57    30 20064.64       52
## attr(,"snr.min")
## [1] 3
## attr(,"scale.range")
## [1]  1 28
## attr(,"length.min")
## [1] 10
## attr(,"noise.span")
## [1] 5
## attr(,"noise.fun")
## [1] "quantile"
## attr(,"noise.min")
## 5% 
## 4.121621

So the main peak close to 50 is correctly found, and the routine picks up another smaller peak at the start.

Upvotes: 1

Related Questions