sadII_lab13_mixedmodels_EM

2241 days ago by agorska

library(flexmix) # http://cran.r-project.org/web/packages/flexmix/vignettes/flexmix-intro.pdf # http://en.wikibooks.org/wiki/Data_Mining_Algorithms_In_R/Clustering/Expectation_Maximization_(EM) 
       
Loading required package: lattice
Loading required package: modeltools
Loading required package: stats4
Loading required package: multcomp
Loading required package: survival
Loading required package: splines
Warning message:
closing unused connection 3
(/srv/nbuser/.sage//temp/sage/20752//interface//tmp20753) 
data("NPreg") head(NPreg) 
       
         x        yn yp yb class id1 id2
1 4.176633 22.380379  4  0     1   1   1
2 1.201631  5.111575  3  0     1   1   1
3 2.295006 13.251058  9  0     1   2   1
4 5.965868 30.285240  3  1     1   2   1
5 2.358083 14.764508  4  0     1   3   2
6 7.637061 42.833760  1  1     1   3   2
png(w=1200) xyplot(yn~x,data=NPreg) xyplot(yp~x,data=NPreg) xyplot(yb~x,data=NPreg) dev.off() 
       
png 
  2 


png(w=1200) xyplot(yn+yp+yb~x,data=NPreg,groups=class, scale=list(y="free")) dev.off() 
       
png 
  2 
?flexmix 
       
WARNING: Output truncated!  
full_output.txt



flexmix                package:flexmix                 R
Documentation

_F_l_e_x_i_b_l_e _M_i_x_t_u_r_e
_M_o_d_e_l_i_n_g

_D_e_s_c_r_i_p_t_i_o_n:

     FlexMix implements a general framework for finite mixtures of
     regression models. Parameter estimation is performed using the
EM
     algorithm: the E-step is implemented by ‘flexmix’, while the
user
     can specify the M-step.

_U_s_a_g_e:

     flexmix(formula, data = list(), k = NULL, cluster = NULL, 
             model=NULL, concomitant=NULL, control = NULL,
             weights = NULL)
     ## S4 method for signature 'flexmix'
     summary(object, eps=1e-4, ...)
     
_A_r_g_u_m_e_n_t_s:

 formula: A symbolic description of the model to be fit. The general
          form is ‘y~x|g’ where ‘y’ is the response, ‘x’ the set of
          predictors and ‘g’ an optional grouping factor for
repeated
          measurements.

    data: An optional data frame containing the variables in the
model.

       k: Number of clusters (not needed if ‘cluster’ is specified).

 cluster: Either a matrix with ‘k’ columns of initial cluster
          membership probabilities for each observation; or a factor
or
          integer vector with the initial cluster assignments of
          observations at the start of the EM algorithm. Default is
          random assignment into ‘k’ clusters.

 weights: An optional vector of replication weights to be used in
the
          fitting process.  Should be ‘NULL’, an integer vector or a
          formula.

   model: Object of class ‘FLXM’ or list of ‘FLXM’ objects. Default
is
          the object returned by calling ‘FLXMRglm()’.

concomitant: Object of class ‘FLXP’. Default is the object returned
by
          calling ‘FLXPconstant’.

 control: Object of class ‘FLXcontrol’ or a named list.

  object: Object of class ‘flexmix’.

     eps: Probabilities below this threshold are treated as zero in
the
          summary method.

     ...: Currently not used.

_D_e_t_a_i_l_s:

     FlexMix models are described by objects of class ‘FLXM’, which
in
     turn are created by driver functions like ‘FLXMRglm’ or

...


     Friedrich Leisch. FlexMix: A general framework for finite
mixture
     models and latent class regression in R. _Journal of
Statistical
     Software_, *11*(8), 2004. http://www.jstatsoft.org/v11/i08/

_S_e_e _A_l_s_o:

     ‘plot-methods’

_E_x_a_m_p_l_e_s:

     data("NPreg", package = "flexmix")
     
     ## mixture of two linear regression models. Note that control
parameters
     ## can be specified as named list and abbreviated if unique.
     ex1 <- flexmix(yn~x+I(x^2), data=NPreg, k=2,
                        control=list(verb=5, iter=100))
     
     ex1
     summary(ex1)
     plot(ex1)
     
     ## now we fit a model with one Gaussian response and one
Poisson
     ## response. Note that the formulas inside the call to FLXMRglm
are
     ## relative to the overall model formula.
     ex2 <- flexmix(yn~x, data=NPreg, k=2,
                    model=list(FLXMRglm(yn~.+I(x^2)), 
                               FLXMRglm(yp~., family="poisson")))
     plot(ex2)
     
     ex2
     table(ex2@cluster, NPreg$class)
     
     ## for Gaussian responses we get coefficients and standard
deviation
     parameters(ex2, component=1, model=1)
     
     ## for Poisson response we get only coefficients
     parameters(ex2, component=1, model=2)
     
     ## fitting a model only to the Poisson response is of course
     ## done like this
     ex3 <- flexmix(yp~x, data=NPreg, k=2,
model=FLXMRglm(family="poisson"))
     
     ## if observations are grouped, i.e., we have several
observations per
     ## individual, fitting is usually much faster:
     ex4 <- flexmix(yp~x|id1, data=NPreg, k=2,
                    model=FLXMRglm(family="poisson"))
     
     ## And now a binomial example. Mixtures of binomials are not
generically
     ## identified, here the grouping variable is necessary:
     set.seed(1234)
     ex5 <- initFlexmix(cbind(yb,1-yb)~x, data=NPreg, k=2,
                        model=FLXMRglm(family="binomial"), nrep=5)
     table(NPreg$class, clusters(ex5))
     
     ex6 <- initFlexmix(cbind(yb,1-yb)~x|id2, data=NPreg, k=2,
                        model=FLXMRglm(family="binomial"), nrep=5)
     table(NPreg$class, clusters(ex6))
     
m1 <- flexmix(yn ~ x + I(x^2), data = NPreg,k=4) #form is ‘y~x|g’ where ‘y’ is the response, ‘x’ the set of # predictors and ‘g’ an optional grouping factor for repeated # measurements. m1 
       
Call:
flexmix(formula = yn ~ x + I(x^2), data = NPreg, k = 4)

Cluster sizes:
 1  2  3  4 
23 74 12 91 

convergence after 40 iterations
summary(m1) 
       
Call:
flexmix(formula = yn ~ x + I(x^2), data = NPreg, k = 4)

       prior size post>0  ratio
Comp.1 0.198   23    141 0.1631
Comp.2 0.297   74    128 0.5781
Comp.3 0.100   12    122 0.0984
Comp.4 0.404   91    130 0.7000

'log Lik.' -635.3649 (df=19)
AIC: 1308.73   BIC: 1371.398 
 
       
names(attributes(m1)) 
       
 [1] "group"       "formula"     "posterior"   "weights"     "iter" 
"cluster"    
 [7] "logLik"      "df"          "control"     "size"       
"converged"   "k0"         
[13] "model"       "prior"       "components"  "concomitant" "call" 
"k"          
[19] "class"      
m1@components 
       
$Comp.1
$Comp.1[[1]]
$coef
(Intercept)           x      I(x^2) 
-2.99380043  6.63776355 -0.09923575 

$sigma
[1] 3.452368



$Comp.2
$Comp.2[[1]]
$coef
(Intercept)           x      I(x^2) 
  1.1082102   3.7332134   0.1174022 

$sigma
[1] 2.482321



$Comp.3
$Comp.3[[1]]
$coef
(Intercept)           x      I(x^2) 
   1.183742   17.820801   -1.779617 

$sigma
[1] 2.520707



$Comp.4
$Comp.4[[1]]
$coef
(Intercept)           x      I(x^2) 
 17.0769218   8.4823893  -0.8371809 

$sigma
[1] 2.677631
xyplot(yn~x,data=NPreg,groups=m1@cluster) xyplot(yn~x,data=NPreg,groups=cut(m1@posterior$scaled[,1],5),auto.key=TRUE) dev.off() 
       
null device 
          1 

# Sprawdzamy statystyczna istotnosc parametrow rm1 <- refit(m1) summary(rm1) 
       
$Comp.1
             Estimate Std. Error z value  Pr(>|z|)    
(Intercept) -2.904157   3.538311 -0.8208 0.4117745    
x            6.568406   1.708574  3.8444 0.0001209 ***
I(x^2)      -0.093885   0.147277 -0.6375 0.5238174    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

$Comp.2
            Estimate Std. Error z value  Pr(>|z|)    
(Intercept) 1.196123   1.315792  0.9091   0.36332    
x           3.673560   0.670818  5.4762 4.345e-08 ***
I(x^2)      0.121803   0.059173  2.0584   0.03955 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

$Comp.3
            Estimate Std. Error  z value Pr(>|z|)    
(Intercept)   1.1539     2.9007   0.3978   0.6908    
x            17.8603     1.4469  12.3436   <2e-16 ***
I(x^2)       -1.7836     0.1401 -12.7308   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

$Comp.4
             Estimate Std. Error z value  Pr(>|z|)    
(Intercept) 17.087153   1.240528  13.774 < 2.2e-16 ***
x            8.476744   0.570653  14.854 < 2.2e-16 ***
I(x^2)      -0.836617   0.053354 -15.681 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
# Teraz inny model rozproszenia - poisson m1p <- flexmix(yp ~ x, data = NPreg, k = 2, model = FLXMRglm(family = "poisson")) summary(m1p) 
       
Call:
flexmix(formula = yp ~ x, data = NPreg, k = 2, model =
FLXMRglm(family = "poisson"))

       prior size post>0 ratio
Comp.1 0.536  112    197 0.569
Comp.2 0.464   88    200 0.440

'log Lik.' -440.6411 (df=5)
AIC: 891.2821   BIC: 907.7737 
lapply(1:2,function(com) parameters(m1p, component = com)) 
       
[[1]]
                     Comp.1
coef.(Intercept)  1.8634875
coef.x           -0.1624198

[[2]]
                     Comp.2
coef.(Intercept) 1.25219374
coef.x           0.06383731
?refit 
       
Generics              package:modeltools               R
Documentation

_G_e_n_e_r_i_c _U_t_i_l_i_t_y
_F_u_n_c_t_i_o_n_s

_D_e_s_c_r_i_p_t_i_o_n:

     A collection of standard generic functions for which other
     packages provide methods.

_U_s_a_g_e:

     ICL(object, ...)
     KLdiv(object, ...)
     Lapply(object, FUN, ...)
     clusters(object, newdata, ...)
     getModel(object, ...)
     parameters(object, ...)
     posterior(object, newdata, ...)
     prior(object, ...)
     refit(object, newdata, ...)
     relabel(object, by, ...)
     
_A_r_g_u_m_e_n_t_s:

  object: S4 classed object.

     FUN: The function to be applied.

 newdata: Optional new data.

      by: Typically a character string specifying how to relabel the
          object.

     ...: Some methods for these generic function may take
additional,
          optional arguments.

_D_e_t_a_i_l_s:

     ICL: Integrated Completed Likelihood criterion for model
          selection.

     KLdiv: Kullback-Leibler divergence.

     Lapply: S4 generic for ‘lapply’

     clusters: Get cluster membership information from a model or
          compute it for new data.

     getModel: Get single model from a collection of models.

     parameters: Get parameters of a model (similar to but more
general
          than ‘coefficients’).

     posterior: Get posterior probabilities from a model or compute
          posteriors for new data.

     prior: Get prior probabilities from a model.

     refit: Refit a model (usually to obtain additional information
          that was not computed or stored during the initial fitting
          process).

     relabel: Relabel a model (usually to obtain a new permutation
of
          labels in mixture models or cluster objects).

_A_u_t_h_o_r(_s):

     Friedrich Leisch
# Narysuj wynik jakosci dopasowania modelu Poissona rm1p <- refit(m1p) summary(rm1p) 
       
$Comp.1
             Estimate Std. Error z value  Pr(>|z|)    
(Intercept)  1.864315   0.174357 10.6925 < 2.2e-16 ***
x           -0.161920   0.039401 -4.1096 3.964e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

$Comp.2
            Estimate Std. Error z value  Pr(>|z|)    
(Intercept) 1.247332   0.179361  6.9543 3.543e-12 ***
x           0.064708   0.027423  2.3596   0.01829 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
############## Grzebiemy w silniku EM # http://cran.r-project.org/web/packages/flexmix/vignettes/flexmix-intro.pdf # str. 13 
       
mymclust <- function ( formula = .~. , diagonal = TRUE ) { require ("mvtnorm") retval <- new ("FLXMC", weighted = TRUE ,formula = formula , dist = "mvnorm", name = "my model - based clustering") retval@defineComponent <- expression ({ #y - dane , logLik <- function (x , y) { #Multivariate Normal distribution density #gęstość rozkładu przy podanych danych dla każdej danej może być inny środek , mean to wektor, # dmvnorm (y , mean = center , sigma = cov , log = TRUE ) } predict <- function (x) { matrix ( center , nrow = nrow (x), ncol = length ( center ), byrow = TRUE ) } new ("FLXcomponent", parameters = list ( center = center , cov = cov ), df = df , logLik = logLik , predict = predict ) }) retval } 
       
?dmvnorm 
       
Mvnorm                 package:mvtnorm                 R
Documentation

_M_u_l_t_i_v_a_r_i_a_t_e _N_o_r_m_a_l
_D_e_n_s_i_t_y _a_n_d _R_a_n_d_o_m
_D_e_v_i_a_t_e_s

_D_e_s_c_r_i_p_t_i_o_n:

     These functions provide the density function and a random
number
     generator for the multivariate normal distribution with mean
equal
     to ‘mean’ and covariance matrix ‘sigma’.

_U_s_a_g_e:

     dmvnorm(x, mean, sigma, log=FALSE)
     rmvnorm(n, mean = rep(0, nrow(sigma)), sigma =
diag(length(mean)),
             method=c("eigen", "svd", "chol"))
     
_A_r_g_u_m_e_n_t_s:

       x: Vector or matrix of quantiles. If ‘x’ is a matrix, each
row
          is taken to be a quantile.

       n: Number of observations.

    mean: Mean vector, default is ‘rep(0, length = ncol(x))’.

   sigma: Covariance matrix, default is ‘diag(ncol(x))’.

     log: Logical; if ‘TRUE’, densities d are given as log(d).

  method: Matrix decomposition used to determine the matrix root of
          ‘sigma’, possible methods are eigenvalue decomposition
          (‘"eigen"’, default), singular value decomposition
(‘"svd"’),
          and Cholesky decomposition (‘"chol"’).

_A_u_t_h_o_r(_s):

     Friedrich Leisch and Fabian Scheipl

_S_e_e _A_l_s_o:

     ‘pmvnorm’, ‘rnorm’, ‘qmvnorm’

_E_x_a_m_p_l_e_s:

     dmvnorm(x=c(0,0))
     dmvnorm(x=c(0,0), mean=c(1,1))
     
     sigma <- matrix(c(4,2,2,3), ncol=2)
     x <- rmvnorm(n=500, mean=c(1,2), sigma=sigma)
     colMeans(x)
     var(x)
     
     x <- rmvnorm(n=500, mean=c(1,2), sigma=sigma, method="chol")
     colMeans(x)
     var(x)
     
     plot(x)
     
data("Nclus") Nclus<-as.data.frame(Nclus) xyplot(V1~V2,data=Nclus) dev.off() 
       
null device 
          1 
m2 <- flexmix(as.matrix(Nclus) ~ 1, k = 4, # Modelujemy dane jako zbior 4ech klastrow o stalych centrach data=Nclus, model = mymclust(diagonal=TRUE)) # Diagonal==TRUE --> model rozproszen jest symetryczny wokol centrow m2 summary(m2) 
       
Warning message:
closing unused connection 4
(/srv/nbuser/.sage//temp/sage/20752//interface//tmp20753) 

Call:
flexmix(formula = as.matrix(Nclus) ~ 1, data = Nclus, k = 4, model =
mymclust(diagonal = TRUE))

Cluster sizes:
  1   2   3   4 
 92  96 213 149 

convergence after 152 iterations

Call:
flexmix(formula = as.matrix(Nclus) ~ 1, data = Nclus, k = 4, model =
mymclust(diagonal = TRUE))

       prior size post>0 ratio
Comp.1 0.159   92    165 0.558
Comp.2 0.175   96    172 0.558
Comp.3 0.397  213    488 0.436
Comp.4 0.269  149    174 0.856

'log Lik.' -2447.3 (df=19)
AIC: 4932.6   BIC: 5014.489 
xyplot(V2~V1,data=Nclus,groups=m2@cluster) dev.off() 
       
null device 
          1 
m3 <- flexmix(as.matrix(Nclus) ~ 1, k = 4, # Modelujemy dane jako zbior 4ech klastrow o stalych centrach data=Nclus, model = mymclust(diagonal=FALSE)) summary(m3) 
       
Call:
flexmix(formula = as.matrix(Nclus) ~ 1, data = Nclus, k = 4, model =
mymclust(diagonal = FALSE))

       prior size post>0 ratio
Comp.1 0.368  203    247 0.822
Comp.2 0.182  100    112 0.893
Comp.3 0.272  150    176 0.852
Comp.4 0.177   97    132 0.735

'log Lik.' -2223.678 (df=23)
AIC: 4493.356   BIC: 4592.484 
xyplot(V2~V1,data=Nclus,groups=m3@cluster) dev.off() 
       
null device 
          1 
# Zadanie: Przerob funkcje mymclust tak aby dopasowywala model z uzyciem zmiennych wyjasniajacych (x) i modelem rozproszenia z rozkladu t-studenta (dmvt) # *: niech df rozkladu t-studenta takze dopasowuje sie podczas ewaluacji modelu. library(mvtnorm) ?dmvt() 
       
# Podpowiedz: do estymacji parametrow chmury wielowymiarowej rozkladu t (zamiast "cov.wt") uzyj funkcji library(sn) mymclust_mle <- function ( formula = .~. ) { require ("sn") require ("mvtnorm") retval <- new ("FLXMC", weighted = TRUE ,formula = formula , dist = "mvtnorm", name = "my t model - based clustering") retval@defineComponent <- expression ({ #y - dane , logLik <- function (x , y) { #Multivariate Normal distribution density #gęstość rozkładu przy podanych danych dla każdej danej może być inny środek , mean to wektor, #dmvt(x, delta, sigma, df = 1, log = TRUE, type = "shifted") dmvt(y , df = df, delta = alpha, sigma= Omega , log = TRUE ) } predict <- function (x) { matrix ( center , nrow = nrow (x), ncol = length ( center ), byrow = TRUE ) } new ("FLXcomponent", #parameters = list ( center = center , cov = cov ), df = df , logLik = logLik , predict = predict ) }) retval@fit <- function (x , y , w ) { b<- mst.mle(y=y, freq = w) para<- b$dp para$center<-para$alpha with (para, eval ( retval@defineComponent ) ) } retval } 
       
m2 <- flexmix(as.matrix(Nclus) ~ 1, k = 4, data=Nclus, model = mymclust_mle()) # Diagonal==TRUE --> model rozproszen jest symetryczny wokol centrow m2 #summary(m2) 
       
Warning messages:
1: closing unused connection 5
(/srv/nbuser/.sage//temp/sage/20752//interface//tmp20753) 
2: closing unused connection 17
(/srv/nbuser/.sage//temp/sage/20752//interface//tmp20753) 
3: closing unused connection 4
(/srv/nbuser/.sage//temp/sage/20752//interface//tmp20753) 
4: closing unused connection 3
(/srv/nbuser/.sage//temp/sage/20752//interface//tmp20753) 


Call:
flexmix(formula = as.matrix(Nclus) ~ 1, data = Nclus, k = 4, model =
mymclust_mle())

Cluster sizes:
  1 
550 

convergence after 11 iterations
summary(m2) 
       
Call:
flexmix(formula = as.matrix(Nclus) ~ 1, data = Nclus, k = 4, model =
mymclust_mle())

       prior size post>0 ratio
Comp.1     1  550    550     1

'log Lik.' -5259.21 (df=16193.17)
AIC: 42904.75   BIC: 112696 
mymclust_mle 
       
Error: object 'mymclust_mle' not found
data(ais, package="sn") attach(ais) X.mat <- model.matrix(~lbm+sex) b <- mst.mle(y=cbind(Ht,Wt)) 
       
The following object(s) are masked from 'ais (position 10)':

    Bfat, bmi, Fe, Hc, Hg, Ht, lbm, rcc, sex, sport, ssf, wcc, Wt
The following object(s) are masked from 'ais (position 12)':

    Bfat, bmi, Fe, Hc, Hg, Ht, lbm, rcc, sex, sport, ssf, wcc, Wt
The following object(s) are masked from 'ais (position 13)':

    Bfat, bmi, Fe, Hc, Hg, Ht, lbm, rcc, sex, sport, ssf, wcc, Wt
The following object(s) are masked from 'ais (position 14)':

    Bfat, bmi, Fe, Hc, Hg, Ht, lbm, rcc, sex, sport, ssf, wcc, Wt
The following object(s) are masked from 'ais (position 15)':

    Bfat, bmi, Fe, Hc, Hg, Ht, lbm, rcc, sex, sport, ssf, wcc, Wt
The following object(s) are masked from 'ais (position 16)':

    Bfat, bmi, Fe, Hc, Hg, Ht, lbm, rcc, sex, sport, ssf, wcc, Wt
The following object(s) are masked from 'ais (position 17)':

    Bfat, bmi, Fe, Hc, Hg, Ht, lbm, rcc, sex, sport, ssf, wcc, Wt
       
$call
mst.mle(y = cbind(Ht, Wt))

$logL
[1] -1452.606

$deviance
[1] 2905.212

$dp
$dp$beta
           Ht       Wt
[1,] 178.7253 64.38539

$dp$Omega
          Ht       Wt
Ht  87.07387 108.3195
Wt 108.31954 271.2117

$dp$alpha
       Ht        Wt 
-2.126663  3.790352 

$dp$df
[1] 18.53857


$se
$se$beta
           Ht       Wt
[1,] 1.690072 2.151787

$se$alpha
       Ht        Wt 
0.4594234 0.7453117 

$se$df
[1] 13.73909

$se$internal
[1] 1.69007161 2.15178718 0.08857694 0.11158461 0.04347682
0.04923446 0.04525679 0.74110818

$se$info
            [,1]       [,2]        [,3]        [,4]         [,5]    
[,6]        [,7]
[1,]   6.5747287 -4.0788936 -21.8810588   0.5955576  -43.1152941  
53.4059075   20.549939
[2,]  -4.0788936  3.7772903   8.3791498  13.2918213   28.6321887  
-0.7652965   31.935287
[3,] -21.8810588  8.3791498 354.8683387 -22.7010305    0.6883544 
-10.3721651   -2.128186
[4,]   0.5955576 13.2918213 -22.7010305 353.9842907   -9.8569654   
6.8292848   13.609755
[5,] -43.1152941 28.6321887   0.6883544  -9.8569654 1109.9214472  
-1.3175621   -8.627341
[6,]  53.4059075 -0.7652965 -10.3721651   6.8292848   -1.3175621
4149.4851608 4119.911358
[7,]  20.5499389 31.9352872  -2.1281863  13.6097551   -8.6273412
4119.9113577 5069.519715
[8,]   0.4491322 -0.3767808 -15.9527991 -16.5938642   -4.7734251   
0.5968082    5.838958
            [,8]
[1,]   0.4491322
[2,]  -0.3767808
[3,] -15.9527991
[4,] -16.5938642
[5,]  -4.7734251
[6,]   0.5968082
[7,]   5.8389578
[8,]   3.8406361


$algorithm
$algorithm$par
[1] 178.7252561  64.3853935   1.8899539   2.8014498  -0.3993911 
-0.2279055   0.2301576   2.9198534

$algorithm$objective
[1] 2905.212

$algorithm$convergence
[1] 0

$algorithm$iterations
[1] 48

$algorithm$evaluations
function gradient 
      65       49 

$algorithm$message
[1] "relative convergence (4)"

$algorithm$value
[1] 2905.212

$algorithm$name
[1] "nlminb"
mymclust <- function ( formula = .~. , diagonal = TRUE ) { require ("mvtnorm") retval <- new ("FLXMC", weighted = TRUE ,formula = formula , dist = "mvnorm", name = "my model - based clustering") retval@defineComponent <- expression ({ #y - dane , logLik <- function (x , y) { #Multivariate Normal distribution density #gęstość rozkładu przy podanych danych dla każdej danej może być inny środek , mean to wektor, #dmvt(x, delta, sigma, df = 1, log = TRUE, type = "shifted") dmvt(y , df = 3, delta = center, sigma = cov , log = TRUE ) } predict <- function (x) { matrix ( center , nrow = nrow (x), ncol = length ( center ), byrow = TRUE ) } new ("FLXcomponent", parameters = list ( center = center , cov = cov ), df = df , logLik = logLik , predict = predict ) }) retval@fit <- function (x , y , w , ...) { para <- cov.wt (y , wt = w )[ c ("center", "cov")] #Weighted covariance df <- (3 * ncol (y) + ncol (y )^2)/2 if ( diagonal ) { para$cov <- diag ( diag ( para$cov )) df <- 2 * ncol (y) } with ( para , eval ( retval@defineComponent )) } retval } 
       
m2 <- flexmix(as.matrix(Nclus) ~ 1, k = 4, data=Nclus, model = mymclust_mle()) # Diagonal==TRUE --> model rozproszen jest symetryczny wokol centrow m2 #summary(m2) 
       
Error in flexmix(as.matrix(Nclus) ~ 1, k = 4, data = Nclus, model =
mymclust_mle()) : 
  error in evaluating the argument 'model' in selecting a method for
function 'flexmix': Error in retval@defineComponent <-
expression({ : 
  object 'retval' not found


Call:
flexmix(formula = as.matrix(Nclus) ~ 1, data = Nclus, k = 4, model =
mymclust(diagonal = TRUE))

Cluster sizes:
  1   2   3   4 
 92  96 213 149 

convergence after 152 iterations