Nothing
R/pFDA-internal.R
In probFDA: Probabilistic Fisher Discriminant Analysis
Defines functions
.pfda.bic
.pfda.fstep
.pfda.loglikelihood
.pfda.main
.pfda.mstep
.pfda.repmat
.pfda.trace
.pfda.bic <-
function(loglik,T,prms,n){
K = prms$K
p = prms$p
comp = switch(prms$model,
'DkBk' = (K-1) + K*p+ (K-1)*(p-K/2) + K^2*(K-1)/2 + K,
'DkB' = (K-1) + K*p+ (K-1)*(p-K/2) + K^2*(K-1)/2 + 1,
'DBk' = (K-1) + K*p+ (K-1)*(p-K/2) + K*(K-1)/2 + K,
'DB' = (K-1) + K*p+ (K-1)*(p-K/2) + K*(K-1)/2 + 1,
'AkjBk'= (K-1) + K*p + (K-1)*(p-K/2) + K^2,
'AkjB' = (K-1) + K*p + (K-1)*(p-K/2) + K*(K-1)+1,
'AkBk' = (K-1) + K*p + (K-1)*(p-K/2) + 2*K,
'AkB' = (K-1) + K*p + (K-1)*(p-K/2) + K+1,
'AjBk' = (K-1) + K*p + (K-1)*(p-K/2) + (K-1)+K,
'AjB' = (K-1) + K*p + (K-1)*(p-K/2) + (K-1)+1,
'ABk' = (K-1) + K*p + (K-1)*(p-K/2) + K+1,
'AB' = (K-1) + K*p + (K-1)*(p-K/2) + 2)
bic = loglik - 1/2 * comp * log(n) # BIC criterion
}
.pfda.fstep <-
function(Y,T,kernel){
n = nrow(Y)
p = ncol(Y)
K = ncol(T)
m = colMeans(Y)
d = min(p-1,(K-1))
# Compute S
XX = as.matrix(Y - t(m*t(matrix(1,n,p))))
TT = t(apply(T,1,"/",sqrt(colSums(T))))
if (n>p & kernel==''){
S = t(XX) %*% XX /n
B = t(XX)%*%TT%*%t(TT)%*%XX / n
#d = rankMatrix(B)
# Eigendecomposition of S^-1 * B
eig = svd(ginv(S)%*%B,nu=d,nv=0) # no sparse case
U = eig$u[,1:d]
}
else{
cat('Kernel mode!\n')
if (n<p | kernel=='linear') G = XX %*% t(XX)
if (kernel=='rbf') {sigma=1; G = as.matrix(exp(dist(XX,diag=T)^2/(2*sigma^2)))}
if (kernel=='sigmoid') {a=1;r=0.1;G = tanh(a * XX %*% t(XX) + r)}
lambda = 0
S = G %*% G + lambda*diag(n)
B = G %*% TT %*% t(TT) %*% G
H = svd(ginv(S)%*%B,nu=d,nv=0)$u[,1:d]
U = svd(t(Y) %*% H,nu=d,nv=0)$u[,1:d]
}
U
}
.pfda.loglikelihood <-
function(prms,Y,V,T){
# Initialization
p = prms$p
K = prms$K
prop = prms$prop
D = prms$D
d = ncol(V)
n = prop * nrow(Y)
# Log-likelihood
ly = 0
for (k in 1:K){ # in the observed space
bk = D[k,p,p]
if (d==1){
Dk = D[k,1,1]
ly = ly - 1/2 * n[k] * (log(Dk) + (p-d)*log(bk) + p*(1+log(2*pi)) - 2*log(prop[k]))
}
else {
Dk = D[k,(1:d),(1:d)] # regularization
ly = ly - 1/2 * n[k] * (log(det(Dk)) + (p-d)*log(bk) + p*(1+log(2*pi)) - 2*log(prop[k]))}
}
ly
}
.pfda.main <-
function(Y,cls,model='AkjBk',kernel='',graph=F){
# Usage: res <- fisherEM(Y,K,...){
# Initialization
n = nrow(Y)
p = ncol(Y)
g = as.numeric(cls)
K = length(unique(g))
T = matrix(0,n,K)
for (i in 1:n){T[i,g[i]] = 1}
#
V = .pfda.fstep(Y,T,kernel)
if (is.matrix(V)) for (i in 1:(K-1)) V[,i] = V[,i] / sqrt(sum(V[,i]^2))
else V = matrix(V / sqrt(sum(V^2)),ncol=1)
prms = .pfda.mstep(Y,V,T,model=model)
Lobs = .pfda.loglikelihood(prms,Y,V,T)
if (graph){
plot(as.data.frame(as.matrix(Y) %*% V[,1:2]),col=max.col(T),xlab='axis 1',ylab='axis 2',pch=20)
}
# Returning the results
crit = .pfda.bic(Lobs,T,prms,n); # BIC
res = list(prms=prms,V=V,bic=crit,ll=Lobs,K=K)
class(res)='pfda'
res
}
.pfda.mstep <-
function(Y,U,T,model){
# 12 different submodels: [DkBk] ... [AkjBk]
# Initialization
Y = as.matrix(Y)
n = nrow(Y)
p = ncol(Y)
K = ncol(T)
d = ncol(U)
mu = matrix(NA,K,K-1)
m = matrix(NA,K,p)
D = array(0,c(K,p,p))
# Projection
X = Y %*% U
nk = colSums(T)
pk = nk/n
# Estimation
for (k in 1:K){
m[k,] = colMeans(Y[T[,k]==1,])
mu[k,] = m[k,] %*% U
YY = as.matrix(Y - .pfda.repmat(m[k,],n,1))
Ck = crossprod(t(.pfda.repmat(T[,k],p,1)) * YY, YY) / (nk[k]-1) #crossprod(x,y) = t(x) %*% y
C = cov(Y)
# Estimation of Delta k amongst 8 submodels
# alpha coefficient:
if (model=='DkBk'){
D[k,(1:d),(1:d)] = crossprod(Ck%*%U,U)
bk = (.pfda.trace(Ck) - sum(diag(crossprod(Ck%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='DkB'){
D[k,(1:d),(1:d)] = crossprod(Ck%*%U,U)
bk = (.pfda.trace(C) - sum(diag(crossprod(C%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='DBk'){
D[k,(1:d),(1:d)] = crossprod(C%*%U,U)
bk = (.pfda.trace(Ck) - sum(diag(crossprod(Ck%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='DB'){
D[k,(1:d),(1:d)] = crossprod(C%*%U,U)
bk = (.pfda.trace(C) - sum(diag(crossprod(C%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='AkjBk'){
if (d==1){D[k,1,1] = diag(crossprod(Ck%*%U,U))} else {
D[k,(1:d),(1:d)] = diag(diag(crossprod(Ck%*%U,U)))}
bk = (.pfda.trace(Ck) - sum(diag(crossprod(Ck%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='AkjB'){
if (d==1){D[k,1,1] = diag(crossprod(Ck%*%U,U))} else {
D[k,(1:d),(1:d)] = diag(diag(crossprod(Ck%*%U,U)))}
bk = (.pfda.trace(C) - sum(diag(crossprod(C%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='AkBk'){
if (d==1){D[k,1,1] = sum(diag(crossprod(Ck%*%U,U)))/d} else{
D[k,(1:d),(1:d)] = diag(rep(sum(diag(crossprod(Ck%*%U,U)))/d,d))}
bk = (.pfda.trace(Ck) - sum(diag(crossprod(Ck%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='AkB'){
if (d==1){D[k,1,1] = sum(diag(crossprod(Ck%*%U,U)))/d} else{
D[k,(1:d),(1:d)] = diag(rep(sum(diag(crossprod(Ck%*%U,U)))/d,d))}
bk = (.pfda.trace(C) - sum(diag(crossprod(C%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='AjBk'){
if (d==1){D[k,1,1] = diag(crossprod(C%*%U,U))} else {
D[k,(1:d),(1:d)] = diag(diag(crossprod(C%*%U,U)))}
bk = (.pfda.trace(Ck) - sum(diag(crossprod(Ck%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='AjB'){
if (d==1){D[k,1,1] = diag(crossprod(C%*%U,U))} else{
D[k,(1:d),(1:d)] = diag(diag(crossprod(C%*%U,U)))}
bk = (.pfda.trace(C) - sum(diag(crossprod(C%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='ABk'){
if (d==1){D[k,1,1] = sum(diag(crossprod(C%*%U,U)))} else {
D[k,(1:d),(1:d)] = diag(rep(sum(diag(crossprod(C%*%U,U)))/d,d))}
bk = (.pfda.trace(Ck) - sum(diag(crossprod(Ck%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='AB'){
if (d==1){D[k,1,1] = sum(diag(crossprod(C%*%U,U)))} else {
D[k,(1:d),(1:d)] = diag(rep(sum(diag(crossprod(C%*%U,U)))/d,d))}
bk = (.pfda.trace(C) - sum(diag(crossprod(C%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
}
prms = list(K=K,p=p,mean=mu,my=m,prop=pk,D=D,model=model)
}
.pfda.repmat <-
function(v,n,p){
if (p==1) M = matrix(1,n,1) %*% v
else M = matrix(rep(v,n),n,(length(v)*p),byrow='T')
M
}
.pfda.trace <-
function(t){
if (length(t)==1) T=t
else T=sum(diag(t))
T
}
.Random.seed <-
c(403L, 10L, -1989895324L, -1355772249L, -392693763L, 1763848692L,
-1594376486L, -1441251291L, 1310297407L, 1492304582L, -1407176880L,
-1963348109L, 1830973841L, 2008196032L, 621954462L, -691194007L,
1510060955L, 196754090L, 445699276L, -200299121L, 1824278981L,
851389820L, -2032527918L, -792235699L, -758303737L, 67142638L,
1374636744L, -1790827061L, -69599191L, -2144475528L, -1199692890L,
-472259135L, 337551955L, 1694738690L, 800268340L, -1411737161L,
1652411821L, -1121270076L, 25307882L, 1155353941L, -678788369L,
735551734L, 300875808L, 1222820323L, 817726081L, -2088817424L,
1798136718L, -1590911751L, -1035180405L, -1851697478L, 1570036156L,
-1282520321L, -1136808363L, -1904887444L, 409004482L, -126090723L,
-2045863273L, -1414188034L, -769060936L, -1992186981L, -1028593159L,
-1503512248L, 14189686L, -308293199L, -1231716381L, 547069106L,
494809348L, 728005191L, 1465761053L, -1024423020L, 1001520954L,
-314530939L, -818587553L, 1261152998L, -61271696L, 881490579L,
1331913521L, 790643168L, -2062125442L, -1861108471L, 1272778555L,
381830922L, -1708951828L, -1944067217L, -2074061659L, 1733412444L,
280727090L, 1007381805L, -1215833113L, -358376690L, -50766616L,
1532344043L, -1891428471L, 1067246360L, -1103153594L, -1020156255L,
158837043L, 520538274L, -665630828L, -1546705385L, -1884990963L,
-1442978012L, -1903910134L, 209825397L, 1387597263L, 762492630L,
2069052032L, 311551171L, -1204813215L, 1678154448L, 269843502L,
1082167769L, 2124467691L, 1917421658L, -1846900196L, -677958881L,
-1746165899L, 873088524L, 711860194L, 460563645L, 256566071L,
180784094L, 40913176L, -1646887813L, 1708335641L, 1400980136L,
620472406L, 309175505L, 285597827L, 1428779794L, -609229148L,
-501427609L, 1582336957L, -94329548L, -522133862L, -400249627L,
934022783L, 1938249734L, -1765958768L, 73655859L, -587336239L,
1465878400L, 1416256094L, -1217543767L, -1165610917L, -1523750806L,
2096056588L, 1640792911L, 1649205253L, -896418116L, 1253316498L,
-1278737267L, -1038294329L, -1498491986L, 350866952L, -2059214325L,
463565545L, -1936065096L, -1873773978L, -482785535L, -1780502125L,
-1964406206L, -1656627724L, 517974007L, -1489003539L, 368899204L,
897198762L, -2042249195L, 39872687L, -1730779466L, -280660768L,
867145763L, 917096257L, -1586451920L, -318949298L, -1506432839L,
-893854517L, -710057606L, 40246396L, -1125181249L, 1475669397L,
609529260L, -2005276798L, 1516957021L, 944317527L, 991396670L,
-1878856840L, 1475840091L, -2142859335L, 871632392L, 2120555702L,
470715377L, -749947357L, -988453262L, 1471265988L, -137414009L,
860271197L, -1603324844L, 2100600442L, -19866171L, -1175740641L,
-310097242L, 956394544L, -1721757997L, -892195599L, -233560032L,
1829162686L, -474191159L, 1376198779L, -1796399542L, 1082460844L,
1962196399L, 837799269L, -1374414052L, 1487978354L, -1125575699L,
-1656597465L, 1051561038L, -171954008L, -623976789L, -395412535L,
802346968L, -651510138L, 2046019425L, 405240563L, 2125079906L,
-1577126956L, 1143305943L, 87981261L, -1208001372L, 630950584L,
-1752042662L, -1163938416L, -1686929028L, 1826133268L, 143438146L,
-1746224384L, -488571972L, 601318224L, 755308834L, -1507263160L,
1728612364L, 1455194204L, 154995922L, -300521472L, 1265507156L,
1143597544L, 2097389562L, -519479600L, -1996950228L, 1905709300L,
-1537828974L, -661799840L, 1382231308L, 1815256800L, -546603678L,
1005351592L, -846722484L, 452667260L, 137476914L, -1284954512L,
1629842756L, 1869126616L, 507068538L, 802354736L, 1204187836L,
-933697708L, -1867460030L, -492342368L, 1628451004L, -1282954800L,
240825698L, 1318509544L, 1543405516L, 1436210940L, -842391566L,
1962577024L, 730536180L, 520131720L, -1563425478L, 369681232L,
190762092L, 1462284436L, -1339270830L, 321996512L, -1682988852L,
931033248L, -639438462L, 454840552L, 359637228L, 1636384828L,
-458153870L, -971454352L, 2093390372L, -1092787656L, 1328091674L,
715307216L, 1517565500L, 735401108L, -632085374L, -966447168L,
-1096190660L, -563681648L, 1605958690L, -265443320L, -1273908084L,
1662412828L, 1654080658L, -1372201856L, -2080506860L, 1512878888L,
-644658118L, -1831813936L, 592137324L, -250890572L, -722395118L,
1587748832L, -1836018740L, -1690389024L, -1752277982L, -1300490968L,
-1363853044L, 604200508L, 1542329586L, -141641488L, -163731644L,
1271070808L, -582534598L, 1983539440L, -139447620L, -641863148L,
831573186L, -1519630752L, 735270268L, -417583920L, -555705310L,
-162199384L, 1265745292L, -693212740L, 582285362L, -746251456L,
-1622301004L, 1800720712L, -151551302L, -1546348784L, 201894124L,
-1394517932L, 1022275346L, -312614496L, 995200972L, 2032547680L,
624702658L, -252666328L, 1323796908L, -2090863812L, 1830736754L,
-75688144L, -91551708L, 574299704L, 1317614298L, 1698867088L,
-1386451716L, 1992037780L, 1596461890L, -301241856L, -853096388L,
-1893275056L, -1955610462L, -865496248L, 1258045964L, 467325660L,
-1812592558L, 1704692864L, 655028052L, -1636524056L, 49830650L,
179224400L, 1014502444L, 1804407156L, -475375982L, 91085408L,
1920161932L, 1169779040L, 102893794L, 1545012392L, -1404543284L,
1206520316L, -420777038L, 1658657776L, -2130245180L, -1770435240L,
-1164363398L, -608396880L, 66060732L, 78114388L, 1642985922L,
1005726496L, -211256388L, -1012936624L, 692429794L, 1316210536L,
-742719668L, 665480956L, -766341390L, 548171008L, 2103698548L,
412069384L, 1828462010L, -673864240L, -818589716L, -432264300L,
1458770130L, -162202528L, -104671924L, -21935328L, -569312638L,
-829005336L, 139420652L, -2099268036L, -247538062L, 287798128L,
1878076708L, 1245362744L, 1507037594L, 1768048080L, 1576822204L,
374916756L, -1582988670L, 1059542464L, 2032193340L, 858493584L,
-1468463326L, 184702472L, 891486988L, 903682972L, 705274898L,
690484096L, 1810967316L, -24029912L, 1166766138L, 951374800L,
144972396L, 2011765812L, -50816878L, -334159520L, -557937332L,
29412832L, -1319472222L, -312915032L, -992637812L, -1944788036L,
1648843890L, -1596524560L, -1568615356L, -2034616232L, 198131770L,
1523232944L, 708751054L, 1075278843L, -363957155L, 1659349898L,
1203608488L, -1578346223L, -9424893L, 1547383556L, 2047561522L,
-67162745L, -62901295L, -1190916778L, 1545642804L, -2146886075L,
-654161265L, -787921912L, 459673638L, -125611821L, -422763979L,
2120540562L, 995411296L, 860482761L, -781873157L, -1505710132L,
1668325978L, 1471799119L, 1248169497L, 2122705966L, 412135996L,
-356266067L, -1150302121L, -840536672L, 1931624190L, -1870733877L,
855018445L, -1336574214L, -1078217416L, 387718881L, -1037962861L,
-1700963916L, -1198437694L, -88098089L, 1460923169L, -518157594L,
779299812L, 1301498645L, 664376127L, -1673546024L, 103428854L,
62983427L, 1709942341L, -1899815326L, 190772688L, 945357561L,
17805355L, 1192886172L, 216133578L, 2042579135L, 1600633353L,
-2123459810L, 251501260L, 1998066333L, -1086652633L, -1857908400L,
-885514258L, -606207845L, 2045831933L, -96890390L, 1583168136L,
2057833521L, -648135005L, 1949316772L, -1763240110L, -1649983833L,
-457047311L, 2005568054L, 937412308L, -373888411L, -1851234193L,
849103528L, 1845237766L, 1969460531L, -302802539L, -955406862L,
-628295104L, 181179305L, 614944283L, 872564716L, -235016198L,
1251644463L, -943856583L, 1629558350L, 520529564L, -1997628403L,
-1677614985L, 292120704L, 711110366L, 378799531L, 2036160173L,
1015819034L, -201485352L, -724971455L, -1159168269L, -1635245804L,
-472550238L, 1623437751L, -2074004095L, -1973419130L, -1240903228L,
-410904331L, -880084641L, 857191352L, 540726934L, 715104419L,
-1421571483L, -547643518L, -1811909392L, -841020263L, -2108054389L,
-1052725892L, 58678186L, -1931944993L, -48721559L, -315310722L,
-925557396L, -826505795L, -169572729L, 1821721200L, 412100878L,
-207653317L, 372106525L, 760393674L, 1732983016L, -1684399151L,
-1406755389L, 403257028L, -943719054L, 1837183047L, -1867025391L,
-2073873002L, 993282292L, 1576256133L, 983444431L, -162094392L,
-1537596570L, -394339949L, 2140705525L, -282572078L, 283309344L,
287636745L, 1080494651L, -840720500L, 464829978L, -1012267377L,
124600793L, 1319220846L, 2055583228L, 1692815213L, 1431353623L,
1217893600L, -1537779010L, 1033132171L, -53819635L, -1582381254L,
1514094456L, -321915615L, 905340371L, 1778926836L, -2027053598L
)
Try the probFDA package in your browser
Any scripts or data that you put into this service are public.
probFDA documentation built on May 1, 2019, 8:48 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .pfda.bic .pfda.fstep .pfda.loglikelihood .pfda.main .pfda.mstep .pfda.repmat .pfda.trace
.pfda.bic <-
function(loglik,T,prms,n){
K = prms$K
p = prms$p
comp = switch(prms$model,
'DkBk' = (K-1) + K*p+ (K-1)*(p-K/2) + K^2*(K-1)/2 + K,
'DkB' = (K-1) + K*p+ (K-1)*(p-K/2) + K^2*(K-1)/2 + 1,
'DBk' = (K-1) + K*p+ (K-1)*(p-K/2) + K*(K-1)/2 + K,
'DB' = (K-1) + K*p+ (K-1)*(p-K/2) + K*(K-1)/2 + 1,
'AkjBk'= (K-1) + K*p + (K-1)*(p-K/2) + K^2,
'AkjB' = (K-1) + K*p + (K-1)*(p-K/2) + K*(K-1)+1,
'AkBk' = (K-1) + K*p + (K-1)*(p-K/2) + 2*K,
'AkB' = (K-1) + K*p + (K-1)*(p-K/2) + K+1,
'AjBk' = (K-1) + K*p + (K-1)*(p-K/2) + (K-1)+K,
'AjB' = (K-1) + K*p + (K-1)*(p-K/2) + (K-1)+1,
'ABk' = (K-1) + K*p + (K-1)*(p-K/2) + K+1,
'AB' = (K-1) + K*p + (K-1)*(p-K/2) + 2)
bic = loglik - 1/2 * comp * log(n) # BIC criterion
}
.pfda.fstep <-
function(Y,T,kernel){
n = nrow(Y)
p = ncol(Y)
K = ncol(T)
m = colMeans(Y)
d = min(p-1,(K-1))
# Compute S
XX = as.matrix(Y - t(m*t(matrix(1,n,p))))
TT = t(apply(T,1,"/",sqrt(colSums(T))))
if (n>p & kernel==''){
S = t(XX) %*% XX /n
B = t(XX)%*%TT%*%t(TT)%*%XX / n
#d = rankMatrix(B)
# Eigendecomposition of S^-1 * B
eig = svd(ginv(S)%*%B,nu=d,nv=0) # no sparse case
U = eig$u[,1:d]
}
else{
cat('Kernel mode!\n')
if (n<p | kernel=='linear') G = XX %*% t(XX)
if (kernel=='rbf') {sigma=1; G = as.matrix(exp(dist(XX,diag=T)^2/(2*sigma^2)))}
if (kernel=='sigmoid') {a=1;r=0.1;G = tanh(a * XX %*% t(XX) + r)}
lambda = 0
S = G %*% G + lambda*diag(n)
B = G %*% TT %*% t(TT) %*% G
H = svd(ginv(S)%*%B,nu=d,nv=0)$u[,1:d]
U = svd(t(Y) %*% H,nu=d,nv=0)$u[,1:d]
}
U
}
.pfda.loglikelihood <-
function(prms,Y,V,T){
# Initialization
p = prms$p
K = prms$K
prop = prms$prop
D = prms$D
d = ncol(V)
n = prop * nrow(Y)
# Log-likelihood
ly = 0
for (k in 1:K){ # in the observed space
bk = D[k,p,p]
if (d==1){
Dk = D[k,1,1]
ly = ly - 1/2 * n[k] * (log(Dk) + (p-d)*log(bk) + p*(1+log(2*pi)) - 2*log(prop[k]))
}
else {
Dk = D[k,(1:d),(1:d)] # regularization
ly = ly - 1/2 * n[k] * (log(det(Dk)) + (p-d)*log(bk) + p*(1+log(2*pi)) - 2*log(prop[k]))}
}
ly
}
.pfda.main <-
function(Y,cls,model='AkjBk',kernel='',graph=F){
# Usage: res <- fisherEM(Y,K,...){
# Initialization
n = nrow(Y)
p = ncol(Y)
g = as.numeric(cls)
K = length(unique(g))
T = matrix(0,n,K)
for (i in 1:n){T[i,g[i]] = 1}
#
V = .pfda.fstep(Y,T,kernel)
if (is.matrix(V)) for (i in 1:(K-1)) V[,i] = V[,i] / sqrt(sum(V[,i]^2))
else V = matrix(V / sqrt(sum(V^2)),ncol=1)
prms = .pfda.mstep(Y,V,T,model=model)
Lobs = .pfda.loglikelihood(prms,Y,V,T)
if (graph){
plot(as.data.frame(as.matrix(Y) %*% V[,1:2]),col=max.col(T),xlab='axis 1',ylab='axis 2',pch=20)
}
# Returning the results
crit = .pfda.bic(Lobs,T,prms,n); # BIC
res = list(prms=prms,V=V,bic=crit,ll=Lobs,K=K)
class(res)='pfda'
res
}
.pfda.mstep <-
function(Y,U,T,model){
# 12 different submodels: [DkBk] ... [AkjBk]
# Initialization
Y = as.matrix(Y)
n = nrow(Y)
p = ncol(Y)
K = ncol(T)
d = ncol(U)
mu = matrix(NA,K,K-1)
m = matrix(NA,K,p)
D = array(0,c(K,p,p))
# Projection
X = Y %*% U
nk = colSums(T)
pk = nk/n
# Estimation
for (k in 1:K){
m[k,] = colMeans(Y[T[,k]==1,])
mu[k,] = m[k,] %*% U
YY = as.matrix(Y - .pfda.repmat(m[k,],n,1))
Ck = crossprod(t(.pfda.repmat(T[,k],p,1)) * YY, YY) / (nk[k]-1) #crossprod(x,y) = t(x) %*% y
C = cov(Y)
# Estimation of Delta k amongst 8 submodels
# alpha coefficient:
if (model=='DkBk'){
D[k,(1:d),(1:d)] = crossprod(Ck%*%U,U)
bk = (.pfda.trace(Ck) - sum(diag(crossprod(Ck%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='DkB'){
D[k,(1:d),(1:d)] = crossprod(Ck%*%U,U)
bk = (.pfda.trace(C) - sum(diag(crossprod(C%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='DBk'){
D[k,(1:d),(1:d)] = crossprod(C%*%U,U)
bk = (.pfda.trace(Ck) - sum(diag(crossprod(Ck%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='DB'){
D[k,(1:d),(1:d)] = crossprod(C%*%U,U)
bk = (.pfda.trace(C) - sum(diag(crossprod(C%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='AkjBk'){
if (d==1){D[k,1,1] = diag(crossprod(Ck%*%U,U))} else {
D[k,(1:d),(1:d)] = diag(diag(crossprod(Ck%*%U,U)))}
bk = (.pfda.trace(Ck) - sum(diag(crossprod(Ck%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='AkjB'){
if (d==1){D[k,1,1] = diag(crossprod(Ck%*%U,U))} else {
D[k,(1:d),(1:d)] = diag(diag(crossprod(Ck%*%U,U)))}
bk = (.pfda.trace(C) - sum(diag(crossprod(C%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='AkBk'){
if (d==1){D[k,1,1] = sum(diag(crossprod(Ck%*%U,U)))/d} else{
D[k,(1:d),(1:d)] = diag(rep(sum(diag(crossprod(Ck%*%U,U)))/d,d))}
bk = (.pfda.trace(Ck) - sum(diag(crossprod(Ck%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='AkB'){
if (d==1){D[k,1,1] = sum(diag(crossprod(Ck%*%U,U)))/d} else{
D[k,(1:d),(1:d)] = diag(rep(sum(diag(crossprod(Ck%*%U,U)))/d,d))}
bk = (.pfda.trace(C) - sum(diag(crossprod(C%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='AjBk'){
if (d==1){D[k,1,1] = diag(crossprod(C%*%U,U))} else {
D[k,(1:d),(1:d)] = diag(diag(crossprod(C%*%U,U)))}
bk = (.pfda.trace(Ck) - sum(diag(crossprod(Ck%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='AjB'){
if (d==1){D[k,1,1] = diag(crossprod(C%*%U,U))} else{
D[k,(1:d),(1:d)] = diag(diag(crossprod(C%*%U,U)))}
bk = (.pfda.trace(C) - sum(diag(crossprod(C%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='ABk'){
if (d==1){D[k,1,1] = sum(diag(crossprod(C%*%U,U)))} else {
D[k,(1:d),(1:d)] = diag(rep(sum(diag(crossprod(C%*%U,U)))/d,d))}
bk = (.pfda.trace(Ck) - sum(diag(crossprod(Ck%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
if (model=='AB'){
if (d==1){D[k,1,1] = sum(diag(crossprod(C%*%U,U)))} else {
D[k,(1:d),(1:d)] = diag(rep(sum(diag(crossprod(C%*%U,U)))/d,d))}
bk = (.pfda.trace(C) - sum(diag(crossprod(C%*%U,U)))) / (p-d)
bk[bk<=0] = 1e-3
D[k,((d+1):p),((d+1):p)] = diag(rep(bk,p-d))
}
}
prms = list(K=K,p=p,mean=mu,my=m,prop=pk,D=D,model=model)
}
.pfda.repmat <-
function(v,n,p){
if (p==1) M = matrix(1,n,1) %*% v
else M = matrix(rep(v,n),n,(length(v)*p),byrow='T')
M
}
.pfda.trace <-
function(t){
if (length(t)==1) T=t
else T=sum(diag(t))
T
}
.Random.seed <-
c(403L, 10L, -1989895324L, -1355772249L, -392693763L, 1763848692L,
-1594376486L, -1441251291L, 1310297407L, 1492304582L, -1407176880L,
-1963348109L, 1830973841L, 2008196032L, 621954462L, -691194007L,
1510060955L, 196754090L, 445699276L, -200299121L, 1824278981L,
851389820L, -2032527918L, -792235699L, -758303737L, 67142638L,
1374636744L, -1790827061L, -69599191L, -2144475528L, -1199692890L,
-472259135L, 337551955L, 1694738690L, 800268340L, -1411737161L,
1652411821L, -1121270076L, 25307882L, 1155353941L, -678788369L,
735551734L, 300875808L, 1222820323L, 817726081L, -2088817424L,
1798136718L, -1590911751L, -1035180405L, -1851697478L, 1570036156L,
-1282520321L, -1136808363L, -1904887444L, 409004482L, -126090723L,
-2045863273L, -1414188034L, -769060936L, -1992186981L, -1028593159L,
-1503512248L, 14189686L, -308293199L, -1231716381L, 547069106L,
494809348L, 728005191L, 1465761053L, -1024423020L, 1001520954L,
-314530939L, -818587553L, 1261152998L, -61271696L, 881490579L,
1331913521L, 790643168L, -2062125442L, -1861108471L, 1272778555L,
381830922L, -1708951828L, -1944067217L, -2074061659L, 1733412444L,
280727090L, 1007381805L, -1215833113L, -358376690L, -50766616L,
1532344043L, -1891428471L, 1067246360L, -1103153594L, -1020156255L,
158837043L, 520538274L, -665630828L, -1546705385L, -1884990963L,
-1442978012L, -1903910134L, 209825397L, 1387597263L, 762492630L,
2069052032L, 311551171L, -1204813215L, 1678154448L, 269843502L,
1082167769L, 2124467691L, 1917421658L, -1846900196L, -677958881L,
-1746165899L, 873088524L, 711860194L, 460563645L, 256566071L,
180784094L, 40913176L, -1646887813L, 1708335641L, 1400980136L,
620472406L, 309175505L, 285597827L, 1428779794L, -609229148L,
-501427609L, 1582336957L, -94329548L, -522133862L, -400249627L,
934022783L, 1938249734L, -1765958768L, 73655859L, -587336239L,
1465878400L, 1416256094L, -1217543767L, -1165610917L, -1523750806L,
2096056588L, 1640792911L, 1649205253L, -896418116L, 1253316498L,
-1278737267L, -1038294329L, -1498491986L, 350866952L, -2059214325L,
463565545L, -1936065096L, -1873773978L, -482785535L, -1780502125L,
-1964406206L, -1656627724L, 517974007L, -1489003539L, 368899204L,
897198762L, -2042249195L, 39872687L, -1730779466L, -280660768L,
867145763L, 917096257L, -1586451920L, -318949298L, -1506432839L,
-893854517L, -710057606L, 40246396L, -1125181249L, 1475669397L,
609529260L, -2005276798L, 1516957021L, 944317527L, 991396670L,
-1878856840L, 1475840091L, -2142859335L, 871632392L, 2120555702L,
470715377L, -749947357L, -988453262L, 1471265988L, -137414009L,
860271197L, -1603324844L, 2100600442L, -19866171L, -1175740641L,
-310097242L, 956394544L, -1721757997L, -892195599L, -233560032L,
1829162686L, -474191159L, 1376198779L, -1796399542L, 1082460844L,
1962196399L, 837799269L, -1374414052L, 1487978354L, -1125575699L,
-1656597465L, 1051561038L, -171954008L, -623976789L, -395412535L,
802346968L, -651510138L, 2046019425L, 405240563L, 2125079906L,
-1577126956L, 1143305943L, 87981261L, -1208001372L, 630950584L,
-1752042662L, -1163938416L, -1686929028L, 1826133268L, 143438146L,
-1746224384L, -488571972L, 601318224L, 755308834L, -1507263160L,
1728612364L, 1455194204L, 154995922L, -300521472L, 1265507156L,
1143597544L, 2097389562L, -519479600L, -1996950228L, 1905709300L,
-1537828974L, -661799840L, 1382231308L, 1815256800L, -546603678L,
1005351592L, -846722484L, 452667260L, 137476914L, -1284954512L,
1629842756L, 1869126616L, 507068538L, 802354736L, 1204187836L,
-933697708L, -1867460030L, -492342368L, 1628451004L, -1282954800L,
240825698L, 1318509544L, 1543405516L, 1436210940L, -842391566L,
1962577024L, 730536180L, 520131720L, -1563425478L, 369681232L,
190762092L, 1462284436L, -1339270830L, 321996512L, -1682988852L,
931033248L, -639438462L, 454840552L, 359637228L, 1636384828L,
-458153870L, -971454352L, 2093390372L, -1092787656L, 1328091674L,
715307216L, 1517565500L, 735401108L, -632085374L, -966447168L,
-1096190660L, -563681648L, 1605958690L, -265443320L, -1273908084L,
1662412828L, 1654080658L, -1372201856L, -2080506860L, 1512878888L,
-644658118L, -1831813936L, 592137324L, -250890572L, -722395118L,
1587748832L, -1836018740L, -1690389024L, -1752277982L, -1300490968L,
-1363853044L, 604200508L, 1542329586L, -141641488L, -163731644L,
1271070808L, -582534598L, 1983539440L, -139447620L, -641863148L,
831573186L, -1519630752L, 735270268L, -417583920L, -555705310L,
-162199384L, 1265745292L, -693212740L, 582285362L, -746251456L,
-1622301004L, 1800720712L, -151551302L, -1546348784L, 201894124L,
-1394517932L, 1022275346L, -312614496L, 995200972L, 2032547680L,
624702658L, -252666328L, 1323796908L, -2090863812L, 1830736754L,
-75688144L, -91551708L, 574299704L, 1317614298L, 1698867088L,
-1386451716L, 1992037780L, 1596461890L, -301241856L, -853096388L,
-1893275056L, -1955610462L, -865496248L, 1258045964L, 467325660L,
-1812592558L, 1704692864L, 655028052L, -1636524056L, 49830650L,
179224400L, 1014502444L, 1804407156L, -475375982L, 91085408L,
1920161932L, 1169779040L, 102893794L, 1545012392L, -1404543284L,
1206520316L, -420777038L, 1658657776L, -2130245180L, -1770435240L,
-1164363398L, -608396880L, 66060732L, 78114388L, 1642985922L,
1005726496L, -211256388L, -1012936624L, 692429794L, 1316210536L,
-742719668L, 665480956L, -766341390L, 548171008L, 2103698548L,
412069384L, 1828462010L, -673864240L, -818589716L, -432264300L,
1458770130L, -162202528L, -104671924L, -21935328L, -569312638L,
-829005336L, 139420652L, -2099268036L, -247538062L, 287798128L,
1878076708L, 1245362744L, 1507037594L, 1768048080L, 1576822204L,
374916756L, -1582988670L, 1059542464L, 2032193340L, 858493584L,
-1468463326L, 184702472L, 891486988L, 903682972L, 705274898L,
690484096L, 1810967316L, -24029912L, 1166766138L, 951374800L,
144972396L, 2011765812L, -50816878L, -334159520L, -557937332L,
29412832L, -1319472222L, -312915032L, -992637812L, -1944788036L,
1648843890L, -1596524560L, -1568615356L, -2034616232L, 198131770L,
1523232944L, 708751054L, 1075278843L, -363957155L, 1659349898L,
1203608488L, -1578346223L, -9424893L, 1547383556L, 2047561522L,
-67162745L, -62901295L, -1190916778L, 1545642804L, -2146886075L,
-654161265L, -787921912L, 459673638L, -125611821L, -422763979L,
2120540562L, 995411296L, 860482761L, -781873157L, -1505710132L,
1668325978L, 1471799119L, 1248169497L, 2122705966L, 412135996L,
-356266067L, -1150302121L, -840536672L, 1931624190L, -1870733877L,
855018445L, -1336574214L, -1078217416L, 387718881L, -1037962861L,
-1700963916L, -1198437694L, -88098089L, 1460923169L, -518157594L,
779299812L, 1301498645L, 664376127L, -1673546024L, 103428854L,
62983427L, 1709942341L, -1899815326L, 190772688L, 945357561L,
17805355L, 1192886172L, 216133578L, 2042579135L, 1600633353L,
-2123459810L, 251501260L, 1998066333L, -1086652633L, -1857908400L,
-885514258L, -606207845L, 2045831933L, -96890390L, 1583168136L,
2057833521L, -648135005L, 1949316772L, -1763240110L, -1649983833L,
-457047311L, 2005568054L, 937412308L, -373888411L, -1851234193L,
849103528L, 1845237766L, 1969460531L, -302802539L, -955406862L,
-628295104L, 181179305L, 614944283L, 872564716L, -235016198L,
1251644463L, -943856583L, 1629558350L, 520529564L, -1997628403L,
-1677614985L, 292120704L, 711110366L, 378799531L, 2036160173L,
1015819034L, -201485352L, -724971455L, -1159168269L, -1635245804L,
-472550238L, 1623437751L, -2074004095L, -1973419130L, -1240903228L,
-410904331L, -880084641L, 857191352L, 540726934L, 715104419L,
-1421571483L, -547643518L, -1811909392L, -841020263L, -2108054389L,
-1052725892L, 58678186L, -1931944993L, -48721559L, -315310722L,
-925557396L, -826505795L, -169572729L, 1821721200L, 412100878L,
-207653317L, 372106525L, 760393674L, 1732983016L, -1684399151L,
-1406755389L, 403257028L, -943719054L, 1837183047L, -1867025391L,
-2073873002L, 993282292L, 1576256133L, 983444431L, -162094392L,
-1537596570L, -394339949L, 2140705525L, -282572078L, 283309344L,
287636745L, 1080494651L, -840720500L, 464829978L, -1012267377L,
124600793L, 1319220846L, 2055583228L, 1692815213L, 1431353623L,
1217893600L, -1537779010L, 1033132171L, -53819635L, -1582381254L,
1514094456L, -321915615L, 905340371L, 1778926836L, -2027053598L
)
Try the probFDA package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.