4
$\begingroup$

I have a large data set composed of several "independent" data frames like this one

     Tiempo UT1 UT2 UT3 UT4 UT5 UT6 UT7 UT8 UT9
3082  616.2   0   4   0   1   1   0   0   0   0
3083  616.4   0   1   0   0   1   0   0   0   0
3084  616.6   0   0   0   0   1   0   0   0   0
3085  616.8   0   1   0   1   2   0   0   0   0
3086  617.0   0   0   0   1   0   0   0   0   0
3087  617.2   0   0   0   0   4   0   0   0   0
3088  617.4   0   1   0   1   0   0   0   0   0
3089  617.6   0   1   0   1   1   0   0   0   0
3090  617.8   0   0   0   1   1   0   0   0   0
3091  618.0   0   1   0   1   0   0   0   0   0
3092  618.2   0   0   0   2   1   0   0   0   0
3093  618.4   1   1   0   0   0   0   0   0   0
3094  618.6   0   0   0  16   4   0   0   6   5
3095  618.8   0  11   0  23   1   0   0   6   4
3096  619.0   1  16   0  23   1   0   0   0   1
3097  619.2   0   6   0   9   0   0   0   1   0
3098  619.4   0   0   1   4   1   0   0   0   1
3099  619.6   0  18   9  16   0   1   1  15   3
3100  619.8   0  21  23  25   0  21   8  22  16
3101  620.0   0  28   3  17   1  10   3   8   2
3102  620.2   0   3   2   7   0  29   0   5   5
3103  620.4   0   2   1   4   0  23   1   2   3
3104  620.6   0   0   0   3   0  19   0   1   1
3105  620.8   0   0   0   3   0  13   0   0   2
3106  621.0   0   0   0   3   1  11   0   1   1
3107  621.2   0   0   0   1   0   8   1   1   1
3108  621.4   0   0   0   2   0  10   0   2   2
3109  621.6   0   0   0   2   0   5   0   3   2
3110  621.8   0   0   0   2   0   7   1   2   1
     UT10 UT11 UT12 UT13 UT14
3082    5    0    0    0    3
3083    1    0    0    0    2
3084    0    0    0    0    1
3085    0    0    0    0    1
3086    0    0    0    0    2
3087    5    0    0    0    2
3088   13    0    0    0    1
3089    2    0    0    0    1
3090    0    0    0    0    0
3091    0    0    0    0    2
3092    1    0    0    0    1
3093    1    1    1    0    0
3094    4    2    0    0    0
3095    6    7    1    1    1
3096   12    1    0    0    3
3097    4    2    0    1    5
3098   10    3    0    0    3
3099   10    1   12    0    1
3100   15   13   21   15    0
3101   11    6    6   14    2
3102    1    3    4    2    0
3103    8    2    2    1    0
3104    1    2    3    1    0
3105    2    4    1    2    0
3106    2    3    1    2    0
3107    0    0    1    2    0
3108    1    0    1    3    0
3109    0    2    1    2    0
3110    0    2    1    4    0

Here's an example of the plot I use to analyse my data (UT4 vs. Tiempo):

My goal is to fit a multi-peak Gaussian of every column UT$_i$ in order to get the parameters for a generic UT and use it for a further statistical analysis. I've found some ideas here using ksmooth fitting multiple peaks to a dataset and extracting individual peak information in R, but the result I got was a unimodal fit of my data. I hope I've made myself clear with the question. Thank you in advance.

Improve this question
$\endgroup$

3 Answers 3

Reset to default
5
$\begingroup$

I had originally intended to post my solution a week ago, on the main stackoverflow.com site, back when this question was first posed over there, however the version which appeared in that forum was placed on hold while I was actually in the midst of entering my solution. Despite my emphatic protests to the moderators, the "on hold" status was allowed to simply languish for a week without any further attention, only to be closed ultimately without comment or explanation. So unfortunately, I suspect that my answer has probably become somewhat stale news by this point, but I'll still post it anyhow, just in case it may still be of some small use to somebody.

To fit two Gaussians, you may use the nls() function, as in the following example:

library(ggplot2)

# Create partial data frame
Tiempo <- seq(616.2, 621.8, 0.2)
UT4 <- c(1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 2, 0, 16, 23, 23, 9, 4, 16, 25,
         17, 7, 4, 3, 3, 3, 1, 2, 2, 2)
df <- data.frame(Tiempo, UT4)

# Fit to a model consisting of a pair of Gaussians.  Note that nls is kind
# of fussy about choosing "reasonable" starting guesses for the parameters.
# It's even fussier if you use the default algorithm (i.e., Gauss-Newton)
fit <- nls(UT4 ~ (C1 * exp(-(Tiempo-mean1)**2/(2 * sigma1**2)) +
                  C2 * exp(-(Tiempo-mean2)**2/(2 * sigma2**2))), data=df,
                  start=list(C1=20, mean1=619, sigma1=1,
                             C2=20, mean2=620, sigma2=1), algorithm="port")  

# Predict the fitted model to a denser grid of x values
dffit <- data.frame(Tiempo=seq(616.2, 621.8, 0.01))
dffit$UT4 <- predict(fit, newdata=dffit)

# Plot the data with the model superimposed
print(ggplot(df, aes(x=Tiempo, y=UT4)) + geom_point() +
             geom_smooth(data=dffit, stat="identity", color="red", size=1.5))

Please beware that nls() can be somewhat fussy about having good starting parameter values. If they aren't at least reasonably close to the correct search neighborhood, you can end up with convergence failures. In particular, the default Gauss-Newton fit algorithm can be especially fussy about good starting values, and may even return a "singular gradient" error (which appears to be distinct, as an error condition, from convergence failure). I chose algorithm="port" in order to mitigate this problem somewhat.

Output is shown below. Also, final tip: if you want to see the optimal final parameter estimates selected by nls(), just type summary(fit) (substituting, if necessary, whatever alternative object variable name you may have chosen in place of fit, to hold the fit information).

Output of example R code

Improve this answer
$\endgroup$
2
  • $\begingroup$ I love you man. Thank you very much for all the hard work it took. You were the only one who understood my real question. I found the code I posted last week but I was having trouble with parameters. I'm a little bit worried about the fact that I have 40 fits to make for each experiment but honestly I think it's the best way to analyse these data. Thank you again! $\endgroup$
    – Matias Andina
    Commented Apr 15, 2014 at 0:26
  • $\begingroup$ You're welcome, and good luck. BTW, if you really have a lot of fits to do, such that it would be prohibitive to select "reasonable" start values for all of them by hand, keep in mind that R does appear to have a tryCatch() function. One alternative strategy for dealing with bad start values might be to just choose like 500 different combinations at random, and keep trying them, in a loop, until one of the combinations leads to a successful fit; i.e., you only exit the loop if the tryCatch() fails to catch an error. $\endgroup$
    – stachyra
    Commented Apr 15, 2014 at 2:31
3
$\begingroup$

I found something interesting searching for the sum of two gaussians and nonlinear models (nls).

The code is large but looks like this.

# Here P is the same as dnorm (probability density function for normal
# distribution), but other functions could be tried here.
P <- function(x, mean, sd)
{
  variance <- sd^2
  exp(-(x-mean)^2/(2*variance)) / sqrt(2*pi*variance)
}
# integrate(P, -1, 1, mean=0, sd=1) is same as integrate(dnorm, -1, 1, mean=0, sd=1)
# integrate(P,-5,5, mean=0, sd=1)   # should be close to 1.0

# Find "peaks" in array.
# R equivalent of Splus peaks() function
# http://finzi.psych.upenn.edu/R/Rhelp02a/archive/33097.html
# (see efg's posting to R-Help on 8 Feb 2007 about problem with ties.)
#
#peaks(c(1,4,4,1,6,1,5,1,1),3)
#[1] FALSE FALSE  TRUE FALSE  TRUE FALSE  TRUE

peaks <- function(series,span=3)
{
  z <- embed(series, span)
  s <- span%/%2
  v <- max.col(z, "first") == 1 + s   # take first if a tie
  result <- c(rep(FALSE,s),v)
  result <- result[1:(length(result)-s)]
  result
}


# First derivative.  Adjust x values to be center of interval.
# Spacing of x-points need not be uniform
Deriv1 <- function(x,y)
{
  y.prime <- diff(y) / diff(x)
  x.prime <- x[-length(x)] + diff(x)/2
  list(x = x.prime,
       y = y.prime)
}

# "Centered" 2nd-derivative. Spacing of x-points assumed to be uniform.
Deriv2 <- function(x,y)
{
  h <- x[2] - x[1]
  Range <- 2:(length(x)-1)  # Drop first and last points
  list(x = x[Range],
       y = (y[Range+1] - 2*y[Range] + y[Range-1]) / h^2)
}


Delta <- 0.01
x <- seq(0.0,10.0, by=Delta)
y1 <- P(x, 3.75, 0.75)
y2 <- P(x, 6.00, 0.50)
y <- y1 + y2

# Approximate area under curve
#0.01*sum(y1)
#0.01*sum(y2)
#0.01*sum(y)

#integrate(P,0,10, mean=3.75, sd=0.75)
#integrate(P,0,10, mean=6.00, sd=0.50)


##### Gaussians
plot(x,y, type="l", lwd=6,
     main="Sum of Two Gaussians",
     xlab="x", ylab="y")
abline(h=0, v=c(3.75,6.00), lty="dotted")
lines(x,y1, col="red")
lines(x,y2, col="blue")

##### 1st Derivative
# According to Abramowitz & Stegun, inflection points are at +/- sigma.
derivative1 <- Deriv1(x,y)

plot(derivative1$x,derivative1$y, type="l",
     main="1st Derivative", xlab="x", ylab="y'")
abline(h=0, v=0, lty="dotted")

peaks.Deriv1   <- peaks( derivative1$y, span=3)
valleys.Deriv1 <- peaks(-derivative1$y, span=3)

points( derivative1$x[peaks.Deriv1], derivative1$y[peaks.Deriv1],
        pch=19, col="red")
points( derivative1$x[valleys.Deriv1], derivative1$y[valleys.Deriv1],
        pch=19, col="blue")

# Approximate location of peak and valley
derivative1$x[peaks.Deriv1]
derivative1$x[valleys.Deriv1]

s.approx <- (derivative1$x[valleys.Deriv1] - derivative1$x[peaks.Deriv1])/2
s.approx

#### 2nd Derivative
derivative2 <- Deriv2(x,y)

plot(derivative2$x,derivative2$y, type="l",
     main="2nd Derivative", xlab="x", ylab="y''")
abline(h=0, v=0, lty="dotted")

peaks.Deriv2   <- peaks( derivative2$y, span=3)
valleys.Deriv2 <- peaks(-derivative2$y, span=3)

points( derivative2$x[peaks.Deriv2], derivative2$y[peaks.Deriv2],
        pch=19, col="red")
points( derivative2$x[valleys.Deriv2], derivative2$y[valleys.Deriv2],
        pch=19, col="blue")

# Approximate location of mean of normal distribution:
derivative2$x[valleys.Deriv2]
derivative2$y[valleys.Deriv2]

mu.approx <-  derivative2$x[valleys.Deriv2]
mu.approx

# Peaks
derivative2$x[peaks.Deriv2]
derivative2$y[peaks.Deriv2]

# Only keep the first four the correpsond to four s.approx values
mu.approx <- mu.approx[1:4]

mtext("U:/efg/lab/R/MixturesOfDistributions/TwoGaussians.R",
      BOTTOM<-1, cex=0.8, adj=0.0, outer=TRUE)

# Fit distribution
# STUPID package:  cannot solve "exact case" without error.
# MUST have noise in data.


fit1 <- nls(y~(a/b)*exp(-(x-c)^2/(2*b^2))+(d/e)*exp(-(x-f)^2/(2*e^2)),
            start=list(a=(1/sqrt(2*pi)) / s.approx[1], b=s.approx[1], c=mu.approx[1],
                       d=(1/sqrt(2*pi)) / s.approx[2], e=s.approx[2], f=mu.approx[2]),
            trace=TRUE)


# add small random noise
set.seed(17)
y <- y + rnorm(length(y), 1E-5, 1E-4)
fit2 <- nls(y~(a/b)*exp(-(x-c)^2/(2*b^2))+(d/e)*exp(-(x-f)^2/(2*e^2)),
            start=list(a=(1/sqrt(2*pi)) / s.approx[1], b=s.approx[1], c=mu.approx[1],
                       d=(1/sqrt(2*pi)) / s.approx[2], e=s.approx[2], f=mu.approx[2]),
            control=nls.control(tol=1E-5, minFactor=1/1024),
            trace=TRUE)
fit2

One of the graphs I get using it is this

enter image description here

I think it's the thing I'm looking for. However, I'll be working with the code in order to write "y" as each column of my data and rewrite parameters (not an easy thing) and tolerance.

Improve this answer
$\endgroup$
2
$\begingroup$

Have you tried R's built in density estimation plot(density(a$UT4)) Should give you an appropriate density estimate. Otherwise you could look into estimating parameters in a normal mixture model (using EM algorithm or something)

=================================

Perhaps what you are trying to do is not 'density estimation' but rather non-linear regression? Perhaps you would like to fit some splines to your data?

Here are some examples of non-linear fits:

plot(a[,6]~a[,2])
lines(supsmu(a[,6],a[,2]),col='magenta',lwd=2)
lines(supsmu(y=a[,6],x=a[,2]),col='cyan',lwd=2)
lines(lowess(y=a[,6],x=a[,2]),col='green',lwd=2)
lines(ksmooth(y=a[,6],x=a[,2]),col='black',lwd=2)
lines(smooth.spline(y=a[,6],a[,2]),col='yellow',lwd=2)
lines(y=gam(a[,6]~s(a[,2]))$fitted.values,x=a[,2],col='brown',lwd=2)

enter image description here

Improve this answer
$\endgroup$
4
  • $\begingroup$ Benjamin, I've tried density and it doesn't do what I want. Moreover, the plot I got its nothing like a fit of the data. I don't know about EM algorithm "or something"... $\endgroup$
    – Matias Andina
    Commented Apr 6, 2014 at 20:07
  • $\begingroup$ what do you mean by "it doesn't do what I want" $\endgroup$
    – bdeonovic
    Commented Apr 6, 2014 at 20:45
  • $\begingroup$ I want to get a smooth version of the yellow fit you proposed with the corresponding parameters. I was thinking that maybe a sum of two gaussians migth do the job but I'm not sure. $\endgroup$
    – Matias Andina
    Commented Apr 8, 2014 at 1:32
  • $\begingroup$ In these plots I am doing regression, so 'sum of two gaussians' doesn't make sense. There is not much data here for these smoothing functions to work with, so as you try to smooth out the curves you end up loosing the two peaks (like the brown and cyan curves) $\endgroup$
    – bdeonovic
    Commented Apr 8, 2014 at 11:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.