################################################################# ### Hemant Ishwaran ### ------------------------------------------------------------- ### ### Normal mixture modeling using the Dirichlet process and Blocked ### Gibbs Sampling. ### ### -------------------------------------------------------------- ### DESCRIPTION: ### Blocked Gibbs (BG) sampler for estimating normal mixture model ### in which the mean & variance are derived from an unknown mixture ### distribution with a finite (but unknown) number of support points ### ### Uses two different approximations for the Dirichlet process: ### (1) Stick Breaking Truncation (2) Weak Limit Finite Dimensional Prior ### Note that (2) is a nice way to put an upper bound on the number ### of mixture components (see Ishwaran et al 2001 and Example 1 below). ### ### Priors are set as in Ishwaran et al (2002) ### -------------------------------------------------------------- ### REQUIRED ARGUMENTS: ### x vector consisting of the observed data ### -------------------------------------------------------------- ### OPTIONAL ARGUMENTS: ### ### A variance for mean (set to 16 times var(x)) ### a shape parameter for variance (default=2) ### a.scale scale parameter for alpha prior (default=2) ### a.shape shape parameter for alpha prior (default=2) ### alpha value for alpha (default action fits alpha ### as a parameter) ### alpha.step step size used in random walk Metropolis ### for alpha (only needed for weak limit ### approximation) ### ### b scale parameter for variance (default=2) ### B variance for hypermean (default=10*A) ### ### cdfOutput generate posterior analysis of cdf's (default F) ### ### Density just do density estimation (default F) ### ### equalVar set all variances to be equal i.e. ### tau[1]=...=tau[N]=tau (default F) ### ### file.out name of file to send graphical output (by default ### no graphics are produced) ### ### mle.flag compute mle; shuts down all verbose details, ### including plots (i.e. overrides file.out) ### ### N number of terms in random measure (default=100) ### n.iter1 number of burn-in iterations ### n.iter2 number of sampled values after burn-in ### ### ### seed set random generator seed ### stickbreak Use a truncated version of the stick-breaker? ### ### tauTrunc use a uniform[0,tau.max] prior for tau variances ### (can only be used when equalVar=F (default setting)) ### tau.max max value when uniform prior used for tau (default= ### var(x)) ### ### verbose.flag shuts down verbose details ### -------------------------------------------------------------- ### RETURNS: ### Bayesian NPMLE using an AIC penalty. ### Plot of estimated density. Graphics for assessing different ### models. (note that graphics are not generated if mle.flag=T) ### -------------------------------------------------------------- ### DETAILS: ### ### (X(i)|Y(i)) ~ N(Y_1(i),Y_2(i)), i = 1,...,n ### (Y(i)|P) ~ P ### P ~ P(N) ### ### where Y(i)=(Y_1(i),Y_2(i)) is a mean-variance ### parameter and P(N) is a random probability measure ### defined by ### ### sum_(k=1,...N) p(k){Z(k) in .}. ### ### Default is to fit the Dirichlet process using a ### truncated stick-breaking construction ### ### p(k)=V(k)(1-V(1))(1-V(2))...(1-V(k-1)) k=1,2,...N-1 ### ### for V(k) i.i.d. Beta(1,alpha) ### ### However, the Dirichlet process can also be approximated ### using the weak limit approximation ### ### (p(1),...,p(N)) ~ Dirichlet(alpha/N,...,alpha/N) ### ### In both cases, the default selection is N=150, ### although the value for N can be specified ### if desired. In the case where the truncation ### approximation is used, truncation diagnostic is ### based on the mean and variance of ### ### U_N = p(N)+p(N+1)+... ### ### Hierarchical model for parameters: ### ### Let Z(k) = (m(k),tau(k)), k=1,...,N ### ### m(k) ~ N(M,A) ### 1/tau(k) ~ gamma(a,b) ### alpha ~ gamma(a.shape,a.scale) ### M ~ N(0,B) ### ### Note: All variances can be set to be equal i.e. ### tau(1)=...=tau(N)=tau by using the equalVar ### option. In this case ### ### 1/tau ~ gamma(a,b) ### ### ### Note: A uniform prior for tau can be invoked by ### using the tauTrunc=T option. In this case ### ### tau ~ uniform[0,tau.max] ### ### ------------------------------------------------------------- ### EXAMPLES: ### ### Example 1 (Finding the Bayesian NPMLE: Two Different Approaches) ### --------- ### Data consists of measurements of Sodium-Lithium Countertransport ### (SLC) on 190 individuals. It is believed that SLC ### measurements are derived from one of two competing genetic ### models: simple dominance, or additive. In the simple ### dominance model there are 2 alleles A1 and A2 and they form ### two phenotypes (A1A1 and A1A2) and (A2A2). In the additive model ### the two alleles form 3 phenotypes (A1A1), (A1A2) and (A2A2). ### Thus the simple dominance model corresponds to a 2-point mixture ### and the additive model to a 3-point mixture. We would like to ### find out which is the true model. ### ### REFERENCES ### ### Roeder, K. (1994) A graphical technique for determining the ### number of components in a mixture of normals. J. of the Amer. ### Statist. Assoc. 89, 487-495. ### ### x=10*c( ### .467,.430,.192,.192,.293,.160,.164,.126,.328,.202,.282,.328,.247,.132, ### .138,.224,.512,.221,.252,.193,.236,.186,.346,.219,.177,.349,.272,.245, ### .213,.197,.229,.245,.210,.281,.175,.273,.439,.471,.451,.237,.313,.136, ### .245,.391,.349,.158,.252,.416,.232,.183,.254,.195,.141,.151,.073,.300, ### .231,.075,.208,.267,.187,.244,.245,.231,.167,.337,.251,.209,.181,.411, ### .191,.288,.280,.119,.394,.443,.423,.534,.393,.273,.149,.225,.159,.170, ### .329,.183,.262,.250,.179,.329,.253,.270,.310,.321,.333,.284,.380,.222, ### .178,.265,.289,.199,.309,.279,.194,.203,.139,.162,.251,.619,.343,.155, ### .340,.332,.412,.218,.304,.261,.206,.231,.182,.267,.198,.191,.258,.179, ### .197,.188,.202,.150,.201,.255,.293,.255,.189,.414,.292,.253,.168,.295, ### .215,.213,.267,.216,.264,.138,.239,.288,.311,.414,.462,.361,.623,.199, ### .215,.321,.273,.259,.206,.376,.228,.155,.186,.097,.176,.174,.386,.393, ### .198,.243,.326,.250,.590,.461,.361,.321,.236,.139,.316,.313,.263,.180, ### .184,.354,.264,.269,.171,.359,.338,.163) ### ### Roeder (1994) argued that the variances are equal in ### this problem. Therefore, we use the equalVar=T option. ### Because of the large number of iterations used (needed to find the ### NPMLE) we turn off verbose details. This suppresses dynamic ### information like number of iterations completed. Here is the ### call to normalMixture. ### ### set.seed(1) ### normalMixture(x,n.iter1=2000,n.iter2=25000,equalVar=T,verbose.flag=F) ### ### The output is: ### ### $mean ### [1] 2.228980 3.835631 5.998110 ### $var ### [1] 0.3311933 0.3311933 0.3311933 ### $prob ### [1] 0.7827038 0.1996390 0.0176562 ### $loglik ### [1] -255.4526 ### ### The mean, var and prob give the bayesian NPMLE. The loglik is ### the penalized log-likelihood value. ### ### The analysis reveals the presence of 3 mixture components for the ### mean (2.2, 3.8 and 5.9) and the value for the variance is ### roughly 0.33. The additive model is the most likely scenario ### here reflecting the evidence of a 3-point mean mixture model. ### ### The preceding analysis was based on infinite dimensional prior. ### Let's consider what happens if we use a prior whose number of ### components N is bounded by the underlying biology. Here we ### set N=3. To call the finite dimensional prior we set the option ### stickbreak=F. ### ### We use the same seed value as before to make the comparison as ### close as possible. Here is the new call: ### ### set.seed(1) ### normalMixture(x,n.iter1=2000,n.iter2=25000,equalVar=T,stickBreak=F,N=3,verbose.flag=F) ### ### The output is: ### ### $mean ### [1] 2.233603 3.806356 5.740366 ### $var ### [1] 0.3359349 0.3359349 0.3359349 ### $prob ### [1] 0.77861661 0.20088694 0.02049645 ### $loglik ### [1] -255.2673 ### ### ### The results are almost identical to our previous values. Only ### one small difference is that the penalized log-likelihood is slightly ### smaller. This reflects a slight improvement in the analysis by ### taking into account the underlying science. In combination with ### our previous analysis the results provide compelling evidence ### of a 3-component mixture model. ### ### ------------------------------------------------------------- ### ### Example 2 (Finding the Bayesian NPMLE: equal and unequal variances) ### --------- ### Velocities in km/sec of 82 galaxies from 6 well-separated ### conic sections of an `unfilled' survey of the Corona Borealis ### region. Multimodality in such surveys is evidence for voids ### and superclusters in the far universe. ### ### REFERENCES ### ### Roeder, K. (1990) Density estimation with confidence sets ### exemplified by superclusters and voids in galaxies. J. Amer. ### Statist. Assoc. 85, 617-624. ### ### x <- c( ### 9172, 9350, 9483, 9558, 9775, 10227, 10406, 16084, ### 16170, 18419, 18552, 18600, 18927, 19052, 19070, 19330, ### 19343, 19349, 19440, 19473, 19529, 19541, 19547, 19663, ### 19846, 19856, 19863, 19914, 19918, 19973, 19989, 20166, ### 20175, 20179, 20196, 20215, 20221, 20415, 20629, 20795, ### 20821, 20846, 20875, 20986, 21137, 21492, 21701, 21814, ### 21921, 21960, 22185, 22209, 22242, 22249, 22314, 22374, ### 22495, 22746, 22747, 22888, 22914, 23206, 23241, 23263, ### 23484, 23538, 23542, 23666, 23706, 23711, 24129, 24285, ### 24289, 24366, 24717, 24990, 25633, 26690, 26995, 32065, ### 32789, 34279) ### ### (A) First we find the NPMLE assuming an equal variance model: ### ### set.seed(1) ### normalMixture(x,n.iter1=2000,n.iter2=25000,equalVar=T,verbose.flag=F) ### ### The ouput yields a 3-component mixture: ### $mean 21396.251 9476.683 32995.854 ### $prob 0.88502649 0.08484509 0.02983183 ### ### ### (B) Next, let's try an unequal variance model: ### ### set.seed(1) ### normalMixture(x,n.iter1=2000,n.iter2=25000,verbose.flag=F) ### ### The ouput once again yields a 3-component mixture model: ### ### $mean 21189.64 9950.36 31221.10 ### $var 5377921 2112348 15805307 ### $prob 0.83211849 0.11342191 0.05314351 ### ### ### (C) Finally, let's try an unequal variance model, using a flat prior ### for the variance: ### ### set.seed(1) ### normalMixture(x,n.iter1=2000,n.iter2=25000,tauTrunc=T,verbose.flag=F,file.out="galaxy") ### ### Now we find a 5 component mixture model: ### ### $mean 9776.947 22825.921 25382.718 19747.474 16159.313 ### $var 341143.91 1650104.91 19507389.04 426082.41 8711.43 ### $prob 0.06071068 0.34265300 0.17407575 0.38632544 0.02460401 ### ### If we check "galaxy.pdf" for the graphical summary, we see that ### the top 25 models oscillate from 3-5 mixtures (page 4). The 3-component ### models look similar to the NPMLE found in (B). ### ### We generally prefer a flat prior for the variance ### ### -------------------------------------------------------------- ### REFERENCES: ### ### Ishwaran and Zarepour (2000). Markov chain Monte Carlo ### in approximate Dirichlet and beta two-parameter process ### hierarchical models. Biometrika 87, pp 371-390. ### ### Ishwaran and James (2001). Gibbs sampling methods for ### stick-breaking priors. JASA 96, pp 161-173. ### ### Ishwaran H., James, L.F. and Sun, J. (2001). ### Bayesian model selection in finite mixtures by marginal density ### decompositions. JASA 96, 1316-1332. ### ### Ishwaran, H. and James, L.F. (2002). ### Approximate Dirichlet process computing in finite normal ### mixtures: smoothing and prior information. J. Comp. Graph. ### Statist. 11, 508-532. ### ### -------------------------------------------------------------- ### BUGS AND SIDE-EFFECTS: ### 1) Graphics need work. ### 2) ???? ### ------------------------------------------------------------ ### Written by: ### ### Hemant Ishwaran ### Department of Quantitative Health ### Sciences ishwaran@bio.ri.ccf.org ### The Cleveland Clinic Foundation phone: 216-444-9932 ### 9500 Euclid Avenue fax: 216-445-2781 ### Cleveland, OH 44195 ### ### http://www.bio.ri.ccf.org/Resume/Pages/Ishwaran ### ------------------------------------------------------------- ### ### THIS PROGRAM SHOULD NOT BE COPIED, USED, MODIFIED, OR ### DISSEMINATED IN ANY WAY WITHOUT SPECIFIC WRITTEN PERMISSION ### FROM THE AUTHOR. PLEASE GIVE CREDIT, WHERE CREDIT IS DUE. ### ################################################################### normalMixture=function(x,N=100,stickbreak=T, file.out=NULL, n.iter1=750,n.iter2=500, alpha=NULL,alpha.step=0.3, a.shape=2,a.scale=2, A=NULL,B=NULL,a=2,b=2, equalVar=F,tauTrunc=F,tau.max=NULL, seed=NULL,verbose.flag=T,mle.flag=F, ...) { ### ------------------------------------------------------------- ### Some useful functions ### # beta-Stacy distribution # log function for big values # log function which zeroes out small values # log of gamma(alpha,beta) density betaStacy=function(alpha,beta,N){ v=c(rbeta(N-1,alpha,beta),1) pr=c(0,cumsum(-log(1-v[-N])),1e300) pr=exp(-pr[-(N+1)])-exp(-pr[-1]) return(list(p=pr/sum(pr),v=v)) } logB=function(a){ifelse(a>=1e300,300/log10(exp(1)),log(a))} logZ=function(p){ifelse(p<=1e-300,0,log(p))} lgd=function(x,a,b){(a-1)*log(x)-x*b} # calculate kseq (C-interface) # calculate mle (C-interface) # calculate predictive density (C-interface) # update for mean,variance (C-interface) # update for tau when using a uniform prior (C-interface) # sample log-gamma for small alpha (C-interface) kseqnorm.C=function(p,x,m,tau,kseq,seed){ n=as.integer(length(x)) N=as.integer(length(p)) storage.mode(p)="double" storage.mode(x)="double" storage.mode(m)="double" storage.mode(tau)="double" storage.mode(kseq)="integer" storage.mode(seed)="integer" r=.C("kseqnorm",p,x,m,tau,kseq=kseq,n,N,seed=seed) list(kseq=r$kseq,seed=r$seed) } mle.C=function(x,p,m,tau){ n=as.integer(length(x)) N=as.integer(length(p)) llk=as.double(0) storage.mode(p)="double" storage.mode(m)="double" storage.mode(tau)="double" .C("map",x,p,m,tau,llk=llk,n,N)$llk } D2.predC=function(p,z,yset,sigmax){ nbreaks=as.integer(length(yset)) N=as.integer(length(p)) area=vector(length=(length(yset)-1)) storage.mode(p)="double" storage.mode(z)="double" storage.mode(yset)="double" storage.mode(sigmax)="double" storage.mode(area)="double" .C("Dpred2",p,z,yset,area=area,sigmax,nbreaks,N)$area } zUpdateNorm.C=function(x,m,tau,M,kseq,kseqn, equalVar,tauTrunc,A,a,b,seed){ n=as.integer(length(x)) N=as.integer(length(m)) storage.mode(x)="double" storage.mode(m)="double" storage.mode(tau)="double" storage.mode(M)="double" storage.mode(kseq)="integer" storage.mode(kseqn)="integer" storage.mode(equalVar)="character" storage.mode(tauTrunc)="character" storage.mode(A)="double" storage.mode(a)="double" storage.mode(b)="double" storage.mode(seed)="integer" r=.C("zUpdateNorm",equalVar,tauTrunc, x,m=m,tau=tau,M,kseq,kseqn,n, N,A,a,b,seed=seed) list(m=r$m,tau=r$tau,seed=r$seed) } gammat.C=function(a,mn,xlow,seed){ n=length(a) gammat=rep(0,n) if (sum(a < -0.5)>0) stop("incorrect shape parameter") pt1=c(1:n)[a>0] pt2=c(1:n)[a==0] pt3=c(1:n)[a==-0.5] if ((n1<-length(pt1))>0){ gammat[pt1]=mn[pt1]/(qgamma(pgamma(xlow*mn[pt1],a[pt1])+ runif(n1)*(1-pgamma(xlow*mn[pt1],a[pt1])),a[pt1])) } if ((n2<-length(pt2))>0){ tiny=1e-6 a2=rep(tiny,n2) gammat[pt2]=mn[pt2]/(qgamma(pgamma(xlow*mn[pt2],a2)+ runif(n2)*(1-pgamma(xlow*mn[pt2],a2)),a2)) } if ((n3<-length(pt3))>0){ tgamma=rep(0,n3) x1=mn[pt3]*xlow xacc=1e-7 storage.mode(n3)="integer" storage.mode(x1)="double" storage.mode(xacc)="double" storage.mode(tgamma)="double" storage.mode(seed)="integer" r=.C("rGammat",n3,x1,xacc,tgamma=tgamma,seed=seed) gammat[pt3]<- mn[pt3]/r$tgamma seed=r$seed } return(list(gammat=gammat,seed=seed)) } l.gammaC=function(a,nsample,seed){ plog=vector(length=nsample) storage.mode(a)="double" storage.mode(plog)="double" storage.mode(nsample)="integer" storage.mode(seed)="integer" r=.C("rlgamma",a,plog=plog,nsample,seed=seed) list(p.log=r$plog,seed=r$seed) } ### ------------------------------------------------------------- ### General initialize ### Rescale x by standard deviation (rescale back at the end) ### cut.pt=0.5 n=length(x) n.rnd=min(100,n.iter2) rnd.pt=round(seq(1,n.iter2,length=n.rnd)) n.breaks=100 n.mle=25 n.sample=0 count=1 x=sort(x) x.org=x sd.xorg=sd(x.org) x=x.org/sd.xorg if (is.null(seed)) seed=round(rnorm(1,sd=1e5)) yset=seq(min(x)-3,max(x)+3,length=n.breaks) D.matrix=matrix(0,ncol=n.breaks-1,nrow=n.rnd) Mean.matrix=Var.matrix=matrix(0,ncol=n,nrow=n.rnd) mle.M=mle.V=mle.P=matrix(NA,ncol=N,nrow=n.mle) mle.llk <- rep(-9e99,n.mle) burnin.flag=T if (mle.flag==T) verbose.flag=F if (verbose.flag==T) cat("starting burn-in iterations ...","\n") ### ------------------------------------------------------------- ### Initialize Bayesian parameters as in Ishwaran et al 2002 ### if (is.null(A)) A=16 if (is.null(B)) B=10*A if (is.null(alpha)){ alpha.flag=T alpha=runif(n=1,min=2,max=5) } else alpha.flag=F kseq=sample(c(1:N),size=n,replace=T,rgamma(N,1)) m=rnorm(N) if (equalVar==T){ tau=rep(rgamma(1,1),N) equalVar="T" tauTrunc="F" } else{ equalVar="F" if (tauTrunc==T){ if (is.null(tau.max)) tau.max=var(x) tau=runif(N)*tau.max tauTrunc="T" } else{ tau=rgamma(N,1) tauTrunc="F" } } M=rnorm(1) #initialization needed for C-interface ### ------------------------------------------------------------- ### GIBBS SAMPLER ### for (i in 1:(n.iter1+n.iter2)){ ### --------------------------------- ### sample (p,z|kseq,tau,M,x,alpha) kseq.uniq=sort(unique(kseq)) kseq.n=rep(0,N) kseq.n[kseq.uniq]=tapply(kseq,kseq,length) n.clus=length(kseq.uniq) if (stickbreak==F){ alpha.cond=rep(alpha,N)/N+kseq.n s.pt=c(1:N)[alpha.cond=cut.pt] #log-gamma's p.log=rep(0,N) p.log[b.pt]=log(rgamma(length(b.pt),alpha.cond[b.pt])) s.p=l.gammaC(alpha.cond[s.pt],length(s.pt),seed) p.log[s.pt]=s.p$p.log;seed=s.p$seed p=exp(p.log)/sum(exp(p.log)) } else{ alpha.cond=rep(1,N)+kseq.n pv=betaStacy(alpha.cond[-N],(alpha+cumsum(kseq.n[N:2]))[(N-1):1],N) p=pv$p;v=pv$v } zUpdateN=zUpdateNorm.C(x,m,tau,M,kseq,kseqn=kseq.n, equalVar,tauTrunc,A,a,b,seed) m=zUpdateN$m tau=zUpdateN$tau seed=zUpdateN$seed if (tauTrunc=="T"){## Use uniform prior for tau? k.pt=(1:N)[kseq.n>0] if (length(k.pt)>0){ mn=tau[k.pt] tau=runif(N)*tau.max r=gammat.C(a=(kseq.n[k.pt]/2-1),mn,xlow=(1/tau.max),seed) tau[k.pt]=r$gammat seed=r$seed tau[tau<1e-10]=1e-10 } else tau=runif(N)*tau.max } ### --------------------------------- ### sample (M|m) v.c=1/(N/A+1/B) M=rnorm(1,mean=(v.c*sum(m)/A),sd=sqrt(v.c)) ### --------------------------------- ### sample (kseq|p,m,tau,x) kseqnorm=kseqnorm.C(p,x,m,tau,kseq,seed) kseq=kseqnorm$kseq seed=kseqnorm$seed Mean=m[kseq] Var=tau[kseq] ### --------------------------------- ### sample alpha if (alpha.flag==T){ if (stickbreak==F){ accept=0 phi=logB(alpha) #transform alpha phi.n=rnorm(1)*alpha.step+phi alpha.n=exp(phi.n) logRatio=((alpha.n-alpha)*sum(p.log)/N- (alpha.n-alpha)*log(sum(exp(p.log))) +log(gamma(alpha.n))-log(gamma(alpha)) +N*(log(gamma(alpha/N))-log(gamma(alpha.n/N))) +(lgd(alpha.n,a.shape,a.scale)-lgd(alpha,a.shape,a.scale)) +(phi.n-phi)) if (runif(1)<=min(1,exp(logRatio))){ alpha=alpha.n accept=1 } } else #use v to work out p alpha=rgamma(1,N+a.shape-1)/(a.scale-sum(logZ(1-v[-N]))) } ### ------------------------------------------------------------- ### send basic info to terminal ### save output values if (i>n.iter1&&verbose.flag==T){ if (burnin.flag==T){ cat("burn-in completed, now outputing values ...","\n") burnin.flag=F } } if (verbose.flag==T){ if (stickbreak==F&&alpha.flag==T&&n.sample>0) cat(paste("iteration",i,", alpha=",signif(alpha,3)), paste("(",signif(100*accept.total/n.sample,3),"% acceptance)",sep=""), paste(" [",n.clus," clusters]",sep=""),"\n") else cat(paste("iteration",i,", alpha=",signif(alpha,3)), paste(" [",n.clus," clusters]",sep=""),"\n") } if (i>n.iter1){ n.sample=n.sample+1 if (rnd.pt[count]==n.sample){ Mean.matrix[count,]=Mean Var.matrix[count,]=Var if (mle.flag==F) D.matrix[count,]=D2.predC(p,z=m,yset,sigmax=tau) count=count+1 } pt.mle=(1:N)[kseq.n>0] n.mm=length(pt.mle) ### BIC and AIC penalization. Comment out as you please. # t.log=mle.C(x,p=p[pt.mle],m=m[pt.mle],tau=tau[pt.mle])-(n.mm-0.5)*log(n) ##bic t.log=mle.C(x,p=p[pt.mle],m=m[pt.mle],tau=tau[pt.mle])-(2*n.mm-1) ##aic if (sum(t.log>mle.llk)>0){ mle.indx=min((1:n.mle)[t.log>mle.llk]) mle.M[mle.indx,]=mle.V[mle.indx,]=mle.P[mle.indx,]=NA mle.M[mle.indx,1:n.mm]=m[kseq.n>0] mle.V[mle.indx,1:n.mm]=tau[kseq.n>0] mle.P[mle.indx,1:n.mm]=p[kseq.n>0] mle.llk[mle.indx]=t.log } } } ### ------------------------------------------------------------- ### Tidy up mle ### order.mle=rev(order(mle.llk)) mle.M=mle.M[order.mle,]*sd.xorg mle.V=mle.V[order.mle,]*(sd.xorg)^2 mle.P=mle.P[order.mle,] mle.llk=mle.llk[order.mle] ### ------------------------------------------------------------- ### Plots ### if (!is.null(file.out) && mle.flag==F){ palette("default") pdf(paste(file.out,".pdf",sep=""),width=8,height=6) matplot(yset[-n.breaks],t(D.matrix), type="l",xaxt="n",xlab="Data",ylab="Density") axis(side=1,at=pretty(yset[-n.breaks]),cex=1.0) plot(range(Mean.matrix*sd.xorg),c(1,n.iter2),type="n", xlab="Posterior Mean", ylab="Iterations") par(cex=0.5) points(x.org,rep(1,n),pch=16,col=2) par(cex=0.2) for (i in 1:n){ jj=jitter(Mean.matrix[,i],amount=0.1) points(jj*sd.xorg,rnd.pt,type="p",pch=16,col=1) lines(jj*sd.xorg,rnd.pt,type="l",lty=3,col=1) } par(cex=1.0) if (length(unique(x.org)0){ wts=mle.P[j,] wts[is.na(wts)]=0 bg.col=round(1000*mle.V[j,]/max(mle.V[j,],na.rm=T)) symbols(mle.M[j,],rep(j,N),add=T,inches=0.1,circles=(0.5+wts),bg=bg.col) } } dev.off() palette("default") } if (verbose.flag==T) cat("finished","\n") ### ------------------------------------------------------------- ### Return top model ### list(mean=mle.M[1,!is.na(mle.M[1,])], var=mle.V[1,!is.na(mle.V[1,])], prob=mle.P[1,!is.na(mle.V[1,])], loglik=mle.llk[1]) }