Thank you Jean

Basically I know that it is possible, however it is quite sensitive to starting values and those 2 peaks do not provide perfect fit. There is probably third peak between those 2 (based on other results).

But if I put third peak to nls I get an error

> fit <- nls(sig ~ d +
+ a1*exp(-0.5*((time-c1)/b1)^2) +
+ a2*exp(-0.5*((time-c2)/b2)^2)+
+ a3*exp(-0.5*((time-c3)/b3)^2),
+ start=list(a1=12.3, b1=18, c1=315, a2=.5, b2=80, c2=390, a3=3, b3=90, c3=480, d=1.5), data=temp)
Error in nls(sig ~ d + a1 * exp(-0.5 * ((time - c1)/b1)^2) + a2 * exp(-0.5 *  :
  singular gradient

There has to be another option which seems ALS function is based on but the source is rather complicated for my math knowledge.

Anyway, thank for your code.

Petr

From: Adams, Jean [mailto:jvadams@usgs.gov]
Sent: Thursday, March 21, 2013 10:58 PM
To: PIKAL Petr
Cc: r-help@r-project.org
Subject: Re: [R] multiple peak fit

Petr,

You could use nonlinear regression to fit two Guassian peaks.  For example,

fit <- nls(sig ~ d + a1*exp(-0.5*((time-c1)/b1)^2) + a2*exp(-0.5*((time-c2)/b2)^2),
          start=list(a1=15, b1=20, c1=315, a2=4, b2=50, c2=465, d=1.5), data=temp)
plot(temp, type="l")
lines(temp$time, predict(fit), lty=2, col="red")

Jean


On Thu, Mar 21, 2013 at 8:17 AM, PIKAL Petr <petr.pikal@precheza.cz<mailto:petr.pikal@precheza.cz>> wrote:
Hi

I went through some extensive search to find suitable method (package, function) to fit multiple peaks. The best I found is ALS package but it requires rather complicated input structure probably resulting from GC-MS experimental data and seems to be an overkill to my problem.

I have basically simple two column data frame with columns time and sig

dput(temp)
structure(list(time = c(33, 34, 37.33333, 44.66667, 56, 70, 85,
99.66667, 114, 127.33333, 140, 151.66667, 163, 174, 185, 196,
206.66667, 217.33333, 228, 238.66667, 249.33333, 260, 271, 282,
292.66667, 303.66667, 314.33333, 325.33333, 336, 346.66667, 357.33333,
368, 379, 389.66667, 400.33333, 411, 421.66667, 432.33333, 443,
454, 465, 476, 487, 497.66667, 508.33333, 519, 530, 541, 551.66667,
562.33333, 573, 584, 595, 606, 616.66667, 627.33333, 638, 649,
659.66667, 670.33333, 681, 692, 703, 713.66667, 724.33333, 735,
746, 757, 768, 779, 789), sig = c(1.558, 1.5549, 1.5619, 1.5614,
1.5618, 1.6044, 1.6161, 1.6287, 1.6432, 1.6925, 1.7273, 1.6932,
1.669, 1.6863, 1.6962, 1.7186, 1.7513, 1.8325, 1.9181, 2.01,
2.1276, 2.2821, 2.5596, 2.9844, 4.1272, 13.421, 15.422, 14.119,
11.491, 8.8799, 6.7774, 5.6223, 4.8775, 4.3628, 4.0517, 3.9146,
3.8704, 3.9162, 4.0372, 4.0948, 4.2054, 4.1221, 3.9145, 3.5724,
3.2108, 2.8311, 2.4605, 2.1985, 1.9685, 1.8158, 1.7487, 1.692,
1.6565, 1.6374, 1.609, 1.5927, 1.5401, 1.5366, 1.5614, 1.5314,
1.4989, 1.5053, 1.4953, 1.4824, 1.4743, 1.4468, 1.4698, 1.4671,
1.4675, 1.4704, 1.4966)), .Names = c("time", "sig"), row.names = c(NA,
-71L), class = "data.frame")

When it is plotted there are clearly 2 peaks visible.

plot(temp, type="l")

Does anybody know how to decompose those data to two (or more) gaussian (or other) peaks?

I can only think about nlme but before I start from scratch I try to ask others.

Best regards
Petr

______________________________________________
R-help@r-project.org<mailto:R-help@r-project.org> mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


	[[alternative HTML version deleted]]

