Скрытая переменная с моделью смеси Гаусса для вменения недостающих данных

В настоящее время я пытаюсь вменять недостающие данные с помощью модели смеси Гаусса. Мой справочный документ отсюда: http://mlg.eng.cam.ac.uk/zoubin/papers/nips93.pdf

В настоящее время я сосредоточен на двумерном наборе данных с двумя гауссовыми компонентами. Это код для определения веса для каждого гауссовского компонента:

myData = faithful[,1:2];    # the data matrix  
for (i in (1:N)) {
        prob1 = pi1*dmvnorm(na.exclude(myData[,1:2]),m1,Sigma1);   # probabilities of sample points under model 1
        prob2 = pi2*dmvnorm(na.exclude(myData[,1:2]),m2,Sigma2);   # same for model 2
        Z<-rbinom(no,1,prob1/(prob1 + prob2 ))    # Z is latent variable as to assign each data point to the particular component 

        pi1<-rbeta(1,sum(Z)+1/2,no-sum(Z)+1/2)
        if (pi1>1/2) {
          pi1<-1-pi1
          Z<-1-Z
        }
      }

Это мой код для определения недостающих значений:

> whichMissXY<-myData[ which(is.na(myData$waiting)),1:2]
> whichMissXY
   eruptions waiting
11     1.833      NA
12     3.917      NA
13     4.200      NA
14     1.750      NA
15     4.700      NA
16     2.167      NA
17     1.750      NA
18     4.800      NA
19     1.600      NA
20     4.250      NA

Мое ограничение заключается в том, как вменять недостающие данные в «ожидающую» переменную на основе конкретного компонента. Этот код - моя первая попытка вменять недостающие данные с использованием условного среднего вменения. Я знаю, это определенно неправильный путь. Результат не будет лгать конкретному компоненту и производить выбросы.

miss.B2 <- which(is.na(myData$waiting))
for (i in miss.B2) {
    myData[i, "waiting"] <- m1[2] + ((rho * sqrt(Sigma1[2,2]/Sigma1[1,1])) * (myData[i, "eruptions"] - m1[1] ) + rnorm(1,0,Sigma1[2,2]))
    #print(miss.B[i,])  
  }

Я был бы признателен, если бы кто-нибудь мог дать какой-либо совет о том, как улучшить метод вменения, который мог бы работать со скрытой / скрытой переменной через модель гауссовой смеси. заранее спасибо


person Jas    schedule 28.11.2016    source источник
comment
Это полностью зависит от ковариационной структуры, которую вы предполагаете для своей модели смеси. Но общий процесс состоит в том, чтобы иметь два шага ЭМ на итерацию.   -  person Alex W    schedule 03.12.2016


Ответы (1)


Это решение для одного типа ковариационной структуры.

devtools::install_github("alexwhitworth/emclustr")
library(emclustr)
data(faithful)
set.seed(23414L)
ff <- apply(faithful, 2, function(j) {
  na_idx <- sample.int(length(j), 50, replace=F)
  j[na_idx] <- NA
  return(j)
})
ff2 <- em_clust_mvn_miss(ff, nclust=2)

# hmm... seems I don't return the imputed values. 
# note to self to update the code    
plot(faithful, col= ff2$mix_est)

введите описание изображения здесь

И выходы параметров

$it
[1] 27

$clust_prop
[1] 0.3955708 0.6044292

$clust_params
$clust_params[[1]]
$clust_params[[1]]$mu
[1]  2.146797 54.833431

$clust_params[[1]]$sigma
[1] 13.41944


$clust_params[[2]]
$clust_params[[2]]$mu
[1]  4.317408 80.398192

$clust_params[[2]]$sigma
[1] 13.71741
person Alex W    schedule 03.12.2016
comment
Дорогой @Alex W, спасибо и жду вашего обновления :) - person Jas; 05.12.2016
comment
@Jas Какое обновление - в моем репо? Для меня это не главный приоритет. Я не планирую в ближайшее время возвращаться к нему. - person Alex W; 05.12.2016