-1

I am trying to fit a non linear function to a given set of data (x and y in code snippet), the function is defined as

f(x) = a/((sin((x-b)/2))^4)

x <- c(0, 5, -5, 10, -10, 15, -15, 20, -20, 25, -25, 30, -30)
y <- c(4.21, 3.73, 2.08, 1.1, 0.61, 0.42, 0.13, 0.1, 0.04, 0.036667, 0.016667, 0.007778, 0.007778)
plot(x,y, log="y")

This is how the initial graph on which I should fit before mentioned function looks like.

Initial graph

But when I try to fit using nls and plot the curve, the graph does not look quite right 

f <- function(x,a,b) { a/((sin((x-b)/2))^4) }
fitmodel <- nls (y ~ f(x,a,b), start=list(a=1,b=1))
lines(x, predict(fitmodel))

This is what I see:

Plotted curve

I am pretty sure I am doing something wrong here and will appreciate any help from you.

alistaire
  • 42,459
  • 4
  • 77
  • 117
Honj25
  • 1
  • 1
  • @ZheyuanLi is right - look at `plot(x,y, log="y", type="l")` for instance to see why this is happening. And then compare `plot(x[order(x)],y[order(x)], log="y", type="l")` – thelatemail Nov 03 '16 at 22:24
  • I think first order y, then x, i.e. `y <- y[order(x)]; x[order(x)]` – Berry Boessenkool Nov 03 '16 at 22:24
  • Inline fix: `lines(sort(x), predict(fitmodel)[order(x)])` – alistaire Nov 03 '16 at 22:32
  • 1
    Thank you for your help, the graph looks definitely better, at least. Unfortunately, the function was given to me, so I could not do much about that. – Honj25 Nov 03 '16 at 22:58

1 Answers1

0

The R interpreter did exactly what you told it to do.

x is unsorted array.

Therefore, predict(fitmodel) make predictions for these unsorted points.

lines(x, predict(fitmodel)) connects the points in the given order. It connects (x[1], predict(fitmodel)[1]) to (x[2], predict(fitmodel)[2]) to (x[3], predict(fitmodel)[3]) etc. Since the points are not sorted by x, you see the picture in the graph.

You might do ind <- order(x); x <- x[ind]; y <- y[ind] as per Zheyuan Li'd suggestion.

Besides, your model makes no sense.

f <- function(x,a,b) { a/((sin((x-b)/2))^4) }
fitmodel <- nls (y ~ f(x,a,b), start=list(a=1,b=1))

For any a and b, f will be a periodic function with period 2π, while your x changes from -30 to 30 with step 5. You cannot reasonably approximate your points with such a function.

user31264
  • 6,557
  • 3
  • 26
  • 40