Nothing
# density function for mixture of normals
dmix <- function(x, alpha,mu1,mu2,sigma1,sigma2) {
if (alpha < 0) return (dnorm(x,mu2,sigma2))
if (alpha > 1) return (dnorm(x,mu1,sigma1))
alpha * dnorm(x,mu1,sigma1) + (1-alpha) * dnorm(x,mu2,sigma2)
}
# log-likelihood
loglik <- function(theta, x) {
alpha <- theta[1]
mu1 <- theta[2]
mu2 <- theta[3]
sigma1 <- theta[4]
sigma2 <- theta[5]
density <- function (x) {
if (alpha < 0) return (Inf)
if (alpha > 1) return (Inf)
if (sigma1<= 0) return (Inf)
if (sigma2<= 0) return (Inf)
dmix(x,alpha,mu1,mu2,sigma1,sigma2)
}
sum( log ( sapply( x, density) ) )
}
# seed the algorithm
m <- mean(faithful$eruptions)
s <- sd(faithful$eruptions)
oldopt <- options(warn=-1) # suppress warnings from log(0)
mle <- nlmax(loglik,p=c(0.5,m-1,m+1,s,s),x=faithful$eruptions)$estimate
mle
options(oldopt)
d <- function(x) {
dmix(x,mle[1],mle[2],mle[3],mle[4],mle[5])
}
myplot <- xhistogram(~eruptions,data=faithful,
n=30,
density=T,
dmath=dmix,
args=list(
alpha=mle[1],
mu1=mle[2],
mu2=mle[3],
sigma1=mle[4],
sigma2=mle[5])
)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.