1

I'm having a problem with my coding. I'm supposed to

  1. Plot a graph of this equation [y=f(x)] where f(x) = (10*((x-1)^2)^(1/3))/(x^2 + 9) for 10001 values of x between (and including) -5 and 5

  2. Among the 10001 x values in a, find the two local maximums of f(x).

I tried doing this : 

# question1
x <- seq(-5,5,length=10001)
y <- (10*((x-1)^2)^(1/3))/(x^2 + 9)

plot(x,y) # This will produce a graph with two max point

# question2
x[which.max(y)]
y[which.max(y)]

however, i only get the coordinates one of the maximum point and clueless as to how do i get another maximum point.

Dave2e
  • 22,192
  • 18
  • 42
  • 50
autumnheartz
  • 11
  • 1
  • 1
  • 4
  • Finding *local maxima* is a common math question. And it is hard to due *well* in a general sense, especially with base R functions. BTW, please be careful to post *actual* R code: your code had three errors (`//` comment, mismatched parens with `{y)`, and `x` used before its definition, as Dave2e was nice enough to find/fix). – r2evans Mar 25 '19 at 16:25
  • `which.max` finds the global maximum, so you can't do it in this way. This may be useful: https://stats.stackexchange.com/questions/22974/how-to-find-local-peaks-valleys-in-a-series-of-data – yarnabrina Mar 25 '19 at 16:26
  • A common (but certainly not infallible) method is to take the derivative (`diff`) and find all tip-over points where the deriv crosses zero (pos-to-neg). (That may be what @yarnabrina's link is doing, I haven't looked at the code closely-enough.) – r2evans Mar 25 '19 at 16:29

1 Answers1

4

You can use find_peaks from package ggpmisc.

library(ggpmisc)
x[ggpmisc:::find_peaks(df$y)]
y[ggpmisc:::find_peaks(df$y)]

The output is:

[1] -1.5  3.0
[1] 1.6373473 0.8818895

Note that the function find_peaks is noted as internal. You therefore need to access it using :::.

You can further parametrize the call to find_peaks using the span and strict argument. See ??find_peaks for details.

You can also directly plot this using the ggplot2 and ggpmisc packages:

x <- seq(-5,5,length=10001)
y <- (10*((x-1)^2)^(1/3))/(x^2 + 9)
df <- data.frame(x = x, y = y)
ggplot(data = df, aes(x = x, y = y)) + geom_line() + stat_peaks(col = "red")

enter image description here

symbolrush
  • 7,123
  • 1
  • 39
  • 67