Questa domanda ha già una risposta qui :

Commenti

  • r non è chiaramente normale. La sua distribuzione è distorta a destra (molti valori piccoli vicini a 0, pochi valori grandi). Lo vedrai digitando " hist (r) ".
  • Suggerimento! Utilizza dput(r) per generare una stringa che può essere facilmente copiata ' n ' incollabile. Ora dobbiamo inserire i dati in r manualmente …
  • @RasmusB å å grazie, io stava cercando quel comando 🙂 Ho modificato la domanda.
  • Non ' tracciare una sequenza di valori di dati per vedere il distribuzione. $ R $ già rappresenta una frequenza o una densità di probabilità dei valori dei dati?
  • Ho risposto a questa domanda più volte in diversi contesti. Una R soluzione per una variabile discreta come Index viene visualizzata in stats.stackexchange. com / a / 43004/919 ; una R soluzione per una variabile continua si trova in stats.stackexchange.com/questions/70153/… ; e una soluzione Excel si trova in stats.stackexchange.com/a/11563/919 .

Risposta

Esiste una differenza tra ladattamento di una distribuzione gaussiana e ladattamento di una curva di densità gaussiana. Cosa normalmixEM sta facendo il primo. Quello che vuoi è (immagino) il secondo.

Adattare una distribuzione è, grosso modo, quello che faresti se facessi un istogramma dei tuoi dati e ho cercato di vedere che tipo di forma avesse. Quello che stai facendo, invece, è semplicemente tracciare una curva. Quella curva ha una gobba nel mezzo, come quello che ottieni tracciando una funzione di densità gaussiana.

Per ottenere ciò che vuoi, devi può utilizzare qualcosa come optim per adattare la curva ai tuoi dati. Il codice seguente utilizzerà i minimi quadrati non lineari per trovare i tre parametri che danno la curva gaussiana più adatta: m è la media gaussiana, s è la deviazione standard e k è un parametro di ridimensionamento arbitrario (poiché la densità gaussiana è vincolata per essere integrata a 1, mentre i tuoi dati non sono “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)) 

Commenti

  • Ho giocato con questa soluzione ma mi hai battuto 🙂 Durante il gioco ho notato che i valori iniziali iniziali dati a optim contavano molto , quindi quando si utilizza questo metodo, assicurarsi di controllare graficamente ladattamento.

Un swer

Propongo di utilizzare minimi quadrati non lineari per questa analisi.

# 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)) 

E dalloutput, sono stato in grado di ottenere la seguente “curva gaussiana” adattata:

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) ) 

inserisci qui la descrizione dellimmagine

Ladattamento non è sorprendente … La funzione $ x \ mapsto \ sin (x) / x $ non sarebbe una modello migliore?

Commenti

  • Grazie. Ottengo la somma dei quadrati residua: 0,01997. Penso di ottenere esattamente lo stesso con la soluzione di Hong Ooi sopra. Lalgoritmo è lo stesso? Inoltre come faccio a tracciare il risultato di nls?
  • Sì, gli algoritmi sono gli stessi, nel senso che se funzionano (don ' t rimanere bloccato in alcuni minimo locale) danno la stessa risposta. A seconda del valore assegnato a method= possono essere esattamente gli stessi.
  • Ho aggiunto due righe per generare il grafico.

Lascia un commento

Il tuo indirizzo email non sarà pubblicato. I campi obbligatori sono contrassegnati *