R/adaptDA-internal.R

Defines functions amdai.crit amdai.estep amdai.initparam amdai.loglikelihood amdai.mstep amdat.crit amdat.estep amdat.estep2 amdat.initiparam amdat.loglikelihood amdat.mstep repmat

.amdai.crit <-
function(loglik,T,prms,n){
  K = prms$K
  C = prms$C
  p = prms$p
  if (prms$model=='qda') comp = (K-C)*p + (K-C)*p*(p+1)/2 + (K-C)
  if (prms$model=='lda') comp = (K-C)*p + p*(p+1)/2 + (K-C)
  aic = 2 * loglik - 2 * comp
  bic = 2 * loglik - comp * log(n)
  icl = 2 * loglik - comp * log(n) + sum(T*log(T))
  crit = list(aic=aic,bic=bic,icl=icl)
}
.amdai.estep <-
function(prms,Y,C,K){
  # Initialization
  ny = nrow(Y)
  p = prms$p
  m = prms$mean
  prop = prms$prop
  V = prms$var
  
  S = matrix(NA,ny,K)
  P = matrix(NA,ny,K)
  
  # Compute posterior probabilities
  for (i in 1:K){
    YY = as.matrix(Y-.repmat(m[i,],ny,1))
    S[,i] = diag(YY %*% solve(V[i,,]) %*% t(YY)) + log(det(V[i,,])) - 2*log(prop[i])
  }
  for (i in 1:K) {P[,i] = 1 / rowSums(exp(0.5*(t(.repmat(S[,i],K,1))-S)))}
  P
}
.amdai.initparam <-
function(obj,Y,K){
  p = ncol(Y)
  ny = nrow(Y)
  C = obj$C
  m = matrix(NA,K,p)
  V = array(NA,c(K,p,p))
  m[1:C,] = obj$mean
  V[1:C,,] = obj$var
  if (K > C){
    for (i in (C+1):K){
      V[i,,] = cov(Y) / K
      m[i,] = mvrnorm(1,colMeans(Y),V[i,,])
    }
  }
  prop = rep(c(1/K),K)
  prms = list(model=obj$model,C=C,K=K,p=p,mean=m,prop=prop,var=V)
}
.amdai.loglikelihood <-
function(prms,T,Y,C,K){
  # Initialization
  ny = nrow(Y)
  p = prms$p
  m = prms$mean
  prop = prms$prop
  V = prms$var
  
  S = matrix(NA,ny,K)
  
  # 	Log-likelihood
  for (i in 1:K){
    YY = as.matrix(Y-.repmat(m[i,],ny,1))
    S[,i] = diag(YY %*% solve(V[i,,]) %*% t(YY)) + log(det(V[i,,])) - 2*log(prop[i]) + p*log(2*pi)
  }
  
  l = - 0.5 * sum(T * S) 
}
.amdai.mstep <-
function(obj,T,Y,C,K){
  # Initialization
  ny = nrow(Y)
  p = ncol(Y)
  m = matrix(NA,K,p)
  prop = rep(c(NA),1,K)
  V = array(NA,c(K,p,p))
  
  # Estimation
  m[1:C,] = obj$mean
  V[1:C,,] = obj$var
  prop[1:C] = obj$prop
  model=obj$model
  for (i in (C+1):K){
    ni = sum(T[,i])
    prop[i] = ni / ny
    m[i,] = colSums(t(.repmat(T[,i],p,1)) * Y) / ni
    XX = as.matrix(Y-.repmat(m[i,],ny,1))
    V[i,,] = t(T[,i] * XX) %*% XX / (ni-1)
  }
  prop[1:C] = (1-sum(prop[(C+1):K])) * prop[1:C]
  
  prms = list(model=model,C=C,K=K,p=p,mean=m,prop=prop,var=V)
}
.amdat.crit <-
function(loglik,T,prms,n){
	K = prms$K
	p = prms$p
	if (prms$model=='qda') comp = K*p + K*p*(p+1)/2 + K -1
	if (prms$model=='lda') comp = K*p + p*(p+1)/2 + K -1
	aic = 2 * loglik - 2 * comp
	bic = 2 * loglik - comp * log(n)
	icl = 2 * loglik - comp * log(n) + sum(T*log(T))
	crit = list(aic=aic,bic=bic,icl=icl)
}
.amdat.estep <-
function(prms,X,S,Y,C,K){
# ne ré-estime pas les t_ik pour X mais uniquement pour X*
	# Initialization
	nx = nrow(X)
	ny = nrow(Y)
	p = prms$p
	m = prms$mean
	prop = prms$prop
	V = prms$V

	T = matrix(NA,ny,K)
	# Compute posterior probabilities of Y
	for (i in 1:K){
		YY = as.matrix(Y-.repmat(m[i,],ny,1))
		T[,i] = diag(YY %*% ginv(V[i,,]) %*% t(YY)) + log(det(V[i,,])) - 2*log(prop[i])
	}
	TT = T
	for (i in 1:K) {TT[,i] = 1 / rowSums(exp(0.5*(t(.repmat(T[,i],K,1))-T)))}
	P = rbind(S,TT)
}
.amdat.estep2 <-
function(prms,X,S,Y,C,K){
# ré-estime les t_ik pour X et X*
	# Initialization
	XY = rbind(X,Y) 
	n = nrow(XY)
	p = ncol(XY)
	m = prms$mean
	prop = prms$prop
	V = prms$V

	T = matrix(NA,n,K)
	# Compute posterior probabilities of Y
	for (i in 1:K){
		YY = as.matrix(XY-.repmat(m[i,],n,1))
		T[,i] = diag(YY %*% ginv(V[i,,]) %*% t(YY)) + log(det(V[i,,])) - 2*log(prop[i])
	}
	P = T
	for (i in 1:K) {P[,i] = 1 / rowSums(exp(0.5*(t(.repmat(T[,i],K,1))-T)))}
	P
}
.amdat.initiparam <-
function(X,cls,Y,K){
	nx = nrow(X)
	p = ncol(X)
	C = max(cls)
	m = matrix(NA,K,p)
	V = array(NA,c(K,p,p))
	for (i in 1:C){
		m[i,] = colMeans(X[cls==i,])
# 		prop[i] = nrow(X[cls==i,]) / n
		V[i,,] = cov(X[cls==i,])
	}
	ny = nrow(Y)
	if (K > C){
		for (i in (C+1):K){
			V[i,,] = cov(Y) / K
			m[i,] = mvrnorm(1,colMeans(Y),V[i,,])
		}
	}
	prop = rep(c(1/K),K)
	prms = list(K=K,p=p,mean=m,prop=prop,V=V)
}
.amdat.loglikelihood <-
function(prms,T,X,Y,C,K){
	# Initialization
	XY = rbind(X,Y) 
	n = nrow(XY)
	p = ncol(XY)
	m = prms$mean
	prop = prms$prop
	V = prms$V

	S = matrix(NA,n,K)

	# Log-likelihood
	for (i in 1:K){
		YY = as.matrix(XY-.repmat(m[i,],n,1))
		S[,i] = diag(YY %*% ginv(V[i,,]) %*% t(YY)) + log(det(V[i,,])) - 2*log(prop[i])
	}

	l = - 0.5 * sum(T * S)
}
.amdat.mstep <-
function(P,X,Y,C,K){
	# Initialization
	p = ncol(X)
	nx = nrow(X)
	ny = nrow(Y)
	n = nx + ny
	m = matrix(NA,K,p)
	prop = rep(c(NA),1,K)
	V = array(NA,c(K,p,p))

	# Estimation
	model = 'qda'
	for (i in 1:K){
		ni = sum(P[1:nx,i])
		nis = sum(P[(nx+1):n,i])
		prop[i] = (ni + nis) / n
		m[i,] = ( colSums(t(.repmat(P[1:nx,i],p,1)) * X) + 
			  colSums(t(.repmat(P[(nx+1):n,i],p,1)) * Y) ) / (ni+nis)
		XX = as.matrix(X-.repmat(m[i,],nx,1))
		XY = as.matrix(Y-.repmat(m[i,],ny,1))
		V[i,,] = (t(P[1:nx,i] * XX) %*% XX + t(P[(nx+1):n,i] * XY) %*% XY) / (ni + nis)
	}
	prms = list(model=model,C=C,K=K,p=p,mean=m,prop=prop,V=V)
}
.Random.seed <-
c(403L, 10L, 1497083353L, -1163633816L, -1834592618L, -1389883375L, 
1544401347L, 1561296338L, 773174372L, 131163047L, -1281864451L, 
1224124660L, -1506574374L, -90746075L, -1263713729L, -1308644922L, 
1767358032L, -1849757069L, 1435928209L, 1811447488L, -1286371170L, 
590750825L, -982904677L, 664146346L, 1756047820L, 848315535L, 
-740795195L, 252775548L, 515326162L, -484663219L, 1372551431L, 
1103654126L, 255903176L, -1961514293L, 134665001L, -140774536L, 
772106406L, -532237631L, 1448330067L, -86224894L, -317539020L, 
1513682615L, 1519159981L, -743289916L, 1244124138L, -1610818475L, 
-2013249041L, 1576569334L, 533810976L, -1760476957L, -2082294399L, 
-1841163280L, -1090102642L, 97545721L, 1219301259L, 1816759226L, 
-1344802116L, -2081887745L, -1699596459L, 1102532716L, 886146242L, 
423685917L, 1472996759L, -1016539906L, 268943032L, -1765793637L, 
-863809799L, 1943388232L, -1564044426L, 155125425L, -1407400221L, 
-1030670926L, 22683140L, -1352854201L, -1006170083L, 1626732180L, 
-609010118L, 467390597L, -313748129L, -1531953690L, -1261480848L, 
-2100590701L, 1188778033L, 212606688L, 1585371006L, 788286985L, 
-1129994693L, 77854730L, 118564844L, 811306095L, 1410032549L, 
1170572124L, -2108575950L, -1300128211L, -319957785L, -1151052274L, 
2067516392L, -1647364629L, -901161335L, 749702680L, -340294330L, 
2066772897L, 1060486707L, 273363874L, 1282721428L, -1135885033L, 
449970957L, -1103121884L, 1615493642L, -2029583499L, 1057998543L, 
-1792095786L, 1822172544L, -1440703037L, 117128033L, 1925620176L, 
-157468370L, -227868455L, 611998443L, 1036994394L, 1730035996L, 
2123188255L, 922621045L, -1952640756L, -332071710L, 1375737277L, 
1844453943L, 1323774174L, -589857768L, -445815941L, 358159129L, 
-44912728L, -882725034L, 1834784721L, -441957501L, 1284220434L, 
-237452892L, 2117891431L, 273382077L, 2068249140L, -1405320806L, 
854875621L, 568473983L, 1684279558L, -158766448L, -320947405L, 
628080849L, 316756608L, 468812638L, -304598871L, 1054430043L, 
-160898710L, -1633254388L, -489429937L, 575853829L, -1927572036L, 
-27451246L, 319771533L, -674356281L, 618846382L, 1120753928L, 
-118200565L, -651009559L, 1254221496L, -1577758362L, 1213764097L, 
93768851L, -1860913342L, 469059828L, 155157751L, -623944467L, 
-352219260L, -1572306518L, -1504430315L, -1401100881L, 392727990L, 
-556141600L, -166880477L, 90348609L, 1238593328L, 1318845774L, 
1541423033L, -1292207669L, 2061728378L, -1482264708L, -782876737L, 
1695066261L, 62477996L, -2071853950L, -2051709859L, 1665952599L, 
1823307326L, -1850066312L, 601278811L, -1257936711L, -1080981240L, 
1868757942L, 42514673L, -2077413597L, 1920842098L, -957256764L, 
-1721137273L, -1472017059L, -572153004L, 733297018L, -1405750075L, 
-561863137L, -513625690L, 1137271600L, 111514067L, 1421531121L, 
1127149856L, -1322250306L, -775946295L, -327381125L, 1100870986L, 
858440108L, -801186641L, -123419035L, 917306908L, -49631118L, 
268822253L, 1416973095L, -1763662002L, 2088308648L, 1662950315L, 
1043339465L, 715989208L, 1091321222L, 766363330L, 156485672L, 
-1810336084L, 440164540L, 2134105586L, -726930128L, 1503904932L, 
-817112392L, -1247192742L, -655308400L, -1257470084L, 18835732L, 
-968450238L, 1986768640L, 125903804L, 1299919696L, 333192994L, 
-1613830840L, -1124061684L, -1522706340L, -1858629422L, -1338847744L, 
1503170900L, -721781272L, -1622705158L, -2118771504L, 1569429804L, 
853436148L, 424017298L, -50569120L, 1506063628L, 546396896L, 
311836002L, 1054925480L, 1120938572L, 123237756L, 1159977266L, 
1068100720L, -613954748L, -894189096L, 1231328378L, 334652976L, 
1181434556L, 150197076L, -388943806L, -656741984L, 1421358780L, 
1982250448L, -1938793118L, 470919144L, -387043380L, -1266359556L, 
140162546L, -535802752L, -2098276108L, 315939464L, -623424710L, 
-306534576L, 1964317292L, -983058284L, 717095250L, 2137659104L, 
871663308L, 333351072L, -1185331326L, -1856854296L, -195809556L, 
-1802825668L, -372693390L, 849654896L, 331762212L, 1572952120L, 
-1328033766L, -586389808L, -269894596L, -976679788L, -1724741502L, 
870821824L, -1927251652L, 1826821776L, 781484066L, -1339475960L, 
1195854988L, -871640036L, -672552302L, -929501056L, 1410159636L, 
-1695069912L, 2118650426L, -247819056L, -2005789588L, 33029812L, 
-460650478L, 1927317472L, 2041648076L, 1787522528L, -294329822L, 
-787940056L, -99281652L, -755790276L, 1305346290L, -766832912L, 
1860423492L, 218177112L, 1390737978L, 1892900592L, 1994371260L, 
637215252L, 1089610946L, 150197344L, 1255120252L, 1497612496L, 
1740010530L, -855650648L, -1133632628L, 1556823484L, 178007090L, 
-1649883328L, 1935135924L, -1230264504L, 1249573562L, 938942224L, 
-558707476L, -17066924L, 327454482L, -1968109664L, -2048904244L, 
1816748896L, 2094077122L, 1450538536L, -892201556L, -89502916L, 
-54134926L, -1233635536L, -2096540124L, -1993141192L, 53326042L, 
-374340208L, -1206457604L, -751909996L, 2059306818L, -532230656L, 
-2079689156L, 1305000528L, 339256994L, -1840643256L, 1613853196L, 
600107228L, 1401662034L, 1417654400L, -964031660L, 1435833832L, 
757391098L, 2016521040L, -649405396L, -1999705740L, -311839086L, 
-302608800L, 103401100L, 1196271456L, -155809054L, -1415343960L, 
-1046713140L, -616422916L, 941690802L, 1278011888L, 1581430212L, 
-1185539240L, -1687982214L, 2064893360L, -972233284L, 1488176212L, 
-1998956606L, -862365920L, -1464816708L, 1798391376L, 1884850658L, 
-119257752L, 1591008588L, -250202884L, -1125651726L, -2104884480L, 
-1684599180L, 517903368L, 74082746L, -610077744L, -293132820L, 
377589652L, 2101870290L, -825836960L, -1983578804L, -1375121120L, 
-217772414L, 522511848L, 1377493996L, -1321067972L, -297220494L, 
-1565459600L, -1790488284L, 78361144L, -871237222L, 33246672L, 
2142932412L, 1689949332L, -1587210110L, 397712832L, 1010503484L, 
-284950384L, 2071694114L, -163100152L, -226907892L, 244631964L, 
-1834038766L, -1715097728L, 552520468L, 1510572328L, -1372338118L, 
-21750832L, -294441364L, -360418764L, 1541006994L, 807820640L, 
545835340L, -1047530528L, 485908898L, -719242840L, 2139453068L, 
-1552945377L, 2035806121L, 499575102L, 1359228972L, -1903903107L, 
-1509154617L, -1553527376L, -906786866L, 930495739L, -807720099L, 
124924042L, 303382184L, 37151249L, 1070764803L, -202529276L, 
1747433010L, -1444256121L, -1389106479L, 436329558L, -1922182092L, 
1681299269L, -1989395057L, -1234801400L, -2057629914L, 94678995L, 
-48185035L, -814057326L, 603453536L, -1179071031L, -1654700805L, 
-573101364L, 2144964442L, -1636230577L, 469671705L, -1734002386L, 
401754940L, -1483312979L, 1985745751L, 896939680L, -2514946L, 
-612927797L, 160266445L, -734503942L, 1634650168L, 1497821665L, 
-1992861037L, -1391095628L, 722796994L, 1786014167L, 2016957985L, 
960515558L, -842660636L, -150242283L, -1709818817L, 1615238104L, 
2054964726L, 2101405187L, -93552315L, -1331878558L, -1351223600L, 
-1622109703L, 990409003L, -2094966116L, 1656500426L, 1736767935L, 
-1223537399L, 1954860574L, -353168436L, 1392492957L, 523074599L, 
-199876016L, 848359150L, 1750149019L, -488565251L, 1170901226L, 
-161498232L, 774059313L, 248343971L, 1029142948L, 194534482L, 
-1663793241L, -1580346383L, -32872650L, -1341340716L, 1263151461L, 
-431917713L, -1729525848L, -397644538L, -475990989L, -1421077867L, 
1836753650L, 1398811456L, 316932265L, 1076938011L, -1482767636L, 
906600186L, 422752047L, -989884615L, 2032111950L, -134336612L, 
-361153267L, -1796987017L, -1626578048L, 345733598L, -192872789L, 
1934721453L, -78211558L, 1930954968L, -1642384063L, 1859108339L, 
-896062444L, 1335685538L, -1458461001L, -851038079L, 2003241094L, 
1086463172L, -1729201675L, -1987251105L, -2132023624L, -693196906L, 
-157762141L, 1953816933L, 1151977602L, -715283984L, -1862633575L, 
-1341542005L, 1102941308L, 840754858L, -842481441L, 1253795433L, 
198046846L, -1333961620L, 1416756413L, -1832249977L, 186773872L, 
-1505532914L, -120088261L, -1896126947L, 266844874L, -62867992L, 
-297592111L, -307730749L, 1919995332L, 1506003570L, -590572729L, 
-310433007L, -2117134698L, 426896372L, -1110993531L, 407750863L, 
1398268360L, 589766758L, -1160284525L, 1614489077L, 404092882L, 
-2065202144L, 1902958089L, -1575923909L, -572325236L, 1305204506L, 
1583774095L, -2097981735L, 1856469358L, -1224588036L, -124232083L, 
821189143L, -971649056L, -123089986L, 707799435L, -1903778941L
)
.repmat <-
function(v,n,p){
	M = matrix(rep(v,n),n,(length(v)*p),byrow=T)
}

Try the adaptDA package in your browser

Any scripts or data that you put into this service are public.

adaptDA documentation built on May 30, 2017, 2:39 a.m.