Existe una diferencia entre ajustar una distribución gaussiana y ajustar una curva de densidad gaussiana. ¿Qué normalmixEM es lo primero. Lo que quieres es (supongo) lo último.
Ajustar una distribución es, en términos generales, lo que harías si hicieras un histograma de sus datos e intentó ver qué tipo de forma tenía. Lo que estás haciendo, en cambio, es simplemente trazar una curva. Esa curva tiene una joroba en el medio, como lo que obtienes al trazar una función de densidad gaussiana.
Para obtener lo que quieres, puede usar algo como optim para ajustar la curva a sus datos. El siguiente código usará mínimos cuadrados no lineales para encontrar los tres parámetros que dan la curva gaussiana que mejor se ajusta: m es la media gaussiana, s es la desviación estándar y k es un parámetro de escala arbitrario (ya que la densidad gaussiana está restringida para integrarse a 1, mientras que sus datos no son «t).
x <- seq_along(r) f <- function(par) { m <- par[1] sd <- par[2] k <- par[3] rhat <- k * exp(-0.5 * ((x - m)/sd)^2) sum((r - rhat)^2) } optim(c(15, 2, 1), f, method="BFGS", control=list(reltol=1e-9))
Comentarios
Propongo utilizar mínimos cuadrados no lineales para este análisis.
# First present the data in a data-frame tab <- data.frame(x=seq_along(r), r=r) #Apply function nls (res <- nls( r ~ k*exp(-1/2*(x-mu)^2/sigma^2), start=c(mu=15,sigma=5,k=1) , data = tab))
Y a partir de la salida, pude obtener la siguiente «curva gaussiana» ajustada:
v <- summary(res)$parameters[,"Estimate"] plot(r~x, data=tab) plot(function(x) v[3]*exp(-1/2*(x-v[1])^2/v[2]^2),col=2,add=T,xlim=range(tab$x) )
El ajuste no es sorprendente … ¿No sería una función «ta $ x \ mapsto \ sin (x) / x $ mejor modelo?
Comentarios