Hello all,
I am trying to fit a truncated mixture model and I wrote a driver for flexmix following the example in the vignette, but it doesn't work for me: it assigns all data points to one component only, e.g.:
> source('bugged.R')

Call:
flexmix(formula = x ~ 1, k = 2, model = truncatedmodel(lower = -4,
    upper = 4))

       prior size post>0 ratio
Comp.1 0.494    0   1000     0
Comp.2 0.506 1000   1000     1

'log Lik.' -707703.3 (df=9)
AIC: 1415425   BIC: 1415469

What am I doing wrong? Please find my code attached.

cheers

--
Giovanni L. Ciampaglia
PhD Student
University of Lugano, MACS Lab

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))
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Reply via email to