library(flexmix)
logf <- .Primitive("log")

dtnorm <- function(
		x,
		lower = -Inf,
		upper = Inf,
		mean = 0,
		sd = 1,
		log = FALSE)
{
	cmin = pnorm(lower,mean,sd)
	cmax = pnorm(upper,mean,sd)
	d = array(0, dim=length(x))
	d[(x >= lower) & (x <= upper)] = dnorm(x,mean,sd) / (cmax - cmin)
	if (log)
		return(logf(d))
	return(d)
}

truncatedmodel <- function (formula = . ~ ., lower = -Inf, upper = Inf) 
{
	require(flexmix)
    z <- new("FLXMC", weighted = TRUE, formula = formula, dist = "tnorm", 
        name = "Truncated bi-modal fitting")
    z@defineComponent <- expression({
        logLik <- function(x, y) 
  			dtnorm(y, lower=lower, upper=upper, mean=mean, sd=sd, log=TRUE)
        predict <- function(x, ...) array(mean, dim = length(x))
        new("FLXcomponent", parameters = list(mean = mean, sd = sd, lower =
			lower, upper = upper), df = df, logLik = logLik, predict = predict)
    })
    z@fit <- function(x, y, w) {
		para <- list( mean = mean(x), sd = sd(x), lower = lower, upper= upper) 
		para$df <- 4
        with(para, eval(z@defineComponent))
    }
    z
}

n <- 1000
x <- c()
while ( length(x) < n) {
	x <- c(x, rnorm( n - length(x) ))
	x <- x[ (x >= -2) & (x <= 2) ]
}
x[1:500] <- x[1:500] - 2
x[501:1000] <- x[501:1000] + 2
m1 <- flexmix( x ~ 1, k = 2, model = truncatedmodel(lower=-4,upper=4))
print(summary(m1))
