linmix_err.pro

basic median filter simulation

                    
                    if ngroup gt 0 then begin
                        
                        xihvar = 1d / (beta[i]^2 / sigsqr[i] + 1d / xixyvar[group] + $
                                       1d / tausqr[k,i])
                        xihat = xihvar * $
                          (xixy[group] / xixyvar[group] + $
                           beta[i] * (eta[xerr[group],i] - alpha[i]) / sigsqr[i] + $
                           mu[k,i] / tausqr[k,i])
                        
                        xi[xerr[group]] = xihat + sqrt(xihvar) * randomn(seed, ngroup)
                        
                    endif
                    
                endfor

            endif
            
            if nyerr gt 0 then begin
                                ;now draw Eta|Xi,x,y,theta
                etaxyvar = yvar[yerr] * (1d - xycorr[yerr]^2)
                etaxy = ygibbs[yerr] + xycov[yerr] / xvar[yerr] * (xi[yerr] - x[yerr])
                if nxnoerr2 gt 0 then etaxy[xnoerr2] = ygibbs[yerr[xnoerr2]]
                etahvar = 1d / (1d / sigsqr[i] + 1d / etaxyvar)
                etahat = etahvar * (etaxy / etaxyvar + $
                                    (alpha[i] + beta[i] * xi[yerr]) / sigsqr[i])
                
                eta[yerr,i] = etahat + sqrt(etahvar) * randomn(seed, nyerr)

            endif

        endif

                                ;now draw new class labels
        if ngauss eq 1 then Glabel[*,i] = 0 else begin
            
            if gibbs then begin
                                ;get unnormalized probability that
                                ;source i came from Gaussian k, given
                                ;xi[i]
                for k = 0, ngauss - 1 do $  
                  gamma[0,k] = pi[k,i] / sqrt(2d * !pi * tausqr[k,i]) * $
                  exp(-0.5 * (xi - mu[k,i])^2 / tausqr[k,i])
                
            endif else begin
                
                for k = 0, ngauss - 1 do begin
                    
                    Sigma11[0,k] = beta[i]^2 * tausqr[k,i] + sigsqr[i] + yvar
                    Sigma12[0,k] = beta[i] * tausqr[k,i] + xycov
                    Sigma22[0,k] = tausqr[k,i] + xvar
                    
                    determ[0, k] = Sigma11[*,k] * Sigma22[*,k] - Sigma12[*,k]^2
                    
                endfor
                
                if ndet gt 0 then begin
                                ;get unnormalized probability that
                                ;source i came from Gaussian k, given
                                ;x[i] and y[i]
                    for k = 0, ngauss - 1 do begin
                        
                        zsqr = (y[det] - alpha[i] - beta[i] * mu[k,i])^2 / Sigma11[det,k] + $
                          (x[det] - mu[k,i])^2 / Sigma22[det,k] - $
                          2d * Sigma12[det,k] * (y[det] - alpha[i] - beta[i] * mu[k,i]) * $
                              (x[det] - mu[k,i]) / (Sigma11[det,k] * Sigma22[det,k])
                        
                        corrz = Sigma12[det,k] / sqrt( Sigma11[det,k] * Sigma22[det,k] )
                        
                        lognorm = -0.5d * alog(determ[det,k]) - 0.5 * zsqr / (1d - corrz^2)
                        
                        gamma[det,k] = pi[k,i] * exp(lognorm) / (2d * !pi)
                            
                    endfor
                    
                endif
                
                if ncens gt 0 then begin
                                ;get unnormalized probability that
                                ;source i came from Gaussian k, given
                                ;x[i] and y[i] > y0[i]
                    for k = 0, ngauss - 1 do begin
                        
                        gamma[cens,k] = pi[k,i] / sqrt(2d * !pi * Sigma22[cens,k]) * $
                          exp(-0.5 * (x[cens] - mu[k,i])^2 / Sigma22[cens,k])
                        
                                ;conditional mean of y, given x
                        cmeany = alpha[i] + beta[i] * mu[k,i] + Sigma12[cens,k] / Sigma22[cens,k] * $
                          (x[cens] - mu[k,i])
                                ;conditional variance of y, given x
                        cvary = Sigma11[cens,k] - Sigma12[cens,k]^2 / Sigma22[cens,k]
                                ;make sure logliky is finite
                        gamma[cens,k] = gamma[cens,k] * gauss_pdf( (y[cens] - cmeany) / sqrt(cvary) )
                        
                    endfor
                    
                endif
                
            endelse
            
            norm = total(gamma, 2)
            
            for j = 0, nx - 1 do begin
                
                gamma0 = reform(gamma[j,*]) / norm[j] ;normalized probability that the i-th data point 
                                                      ;is from the k-th Gaussian, given the observed
                                                      ;data point 
                Gjk = multinom(1, gamma0, seed=seed)
                Glabel[j,i] = where(Gjk eq 1)
                
            endfor
            
        endelse
        
;now draw new values of regression parameters, theta = (alpha, beta,
;sigsqr)
        
        if gibbs then begin
                                ;use gibbs sampler to draw alpha,beta|Xi,Eta,sigsqr
            Xmat = [[replicate(1d, nx)], [xi]]
            Vcoef = invert( Xmat ## transpose(Xmat), /double ) * sigsqr[i]
            
            coefhat = linfit( xi, eta[*,i] )
            coef = coefhat + mrandomn(seed, Vcoef)
            
            alpha[i] = coef[0]
            beta[i] = coef[1]
            
        endif else begin

            theta = [alpha[i], beta[i], sigsqr[i]]

            loglik = loglik_mixerr( x, ygibbs, xvar, yvar, xycov, delta, theta, $
                                    pi[*,i], mu[*,i], tausqr[*,i], Glabel[*,i] )
            logprior = logprior_mixerr(mu[*,i], mu0[i], tausqr[*,i], usqr[i], wsqr[i])

            logpost[i] = loglik + logprior ;log-posterior for current parameter values

                                ;use metropolis update to get new
                                ;values of the coefficients
            coef = [alpha[i], beta[i]] + mrandomn(seed, jvar_coef)
            
            theta = [coef[0], coef[1], sigsqr[i]]
            loglik_new = loglik_mixerr( x, ygibbs, xvar, yvar, xycov, delta, theta, $
                                        pi[*,i], mu[*,i], tausqr[*,i], Glabel[*,i] )
            logprior_new = logprior_mixerr(mu[*,i], mu0[i], tausqr[*,i], usqr[i], wsqr[i])
            
            logpost_new = loglik_new + logprior_new
            
            accept = linmix_metro_update( logpost_new, logpost[i], seed )
            
            if accept then begin
                
                naccept[0] = naccept[0] + 1L
                alpha[i] = coef[0]
                beta[i] = coef[1]
                logpost[i] = logpost_new
                
            endif
            
        endelse
                                ;now get sigsqr
        if gibbs then begin
            
            ssqr = total( (eta[*,i] - alpha[i] - beta[i] * xi)^2 ) / (nx - 2)
            sigsqr[i] = (nx - 2) * ssqr / randomchi(seed, nx - 2.0)
            
        endif else begin
                                ;do metropolis update
            log_ssqr = alog(sigsqr[i]) + sqrt(jvar_ssqr) * randomn(seed)
            ssqr = exp(log_ssqr)
            
            theta = [alpha[i], beta[i], ssqr]
            
            loglik_new = loglik_mixerr( x, ygibbs, xvar, yvar, xycov, delta, theta, $
                                        pi[*,i], mu[*,i], tausqr[*,i], Glabel[*,i] )
            logprior_new = logprior_mixerr(mu[*,i], mu0[i], tausqr[*,i], usqr[i], wsqr[i])
            
            logpost_new = loglik_new + logprior_new + log_ssqr
            logpost_old = logpost[i] + alog(sigsqr[i])
            
            accept = linmix_metro_update( logpost_new, logpost_old, seed )
            
            if accept then begin
                
                naccept[1] = naccept[1] + 1L
                sigsqr[i] = ssqr
                logpost[i] = loglik_new + logprior_new
                
            endif
            
        endelse 

;now do mixture model parameters, psi = (pi,mu,tausqr)
        
        if gibbs then begin

            for k = 0, ngauss - 1 do begin 

                group = where(Glabel[*,i] eq k, ngroup)
                nk[k] = ngroup

                if ngroup gt 0 then begin

                                ;get mu|Xi,G,tausqr,mu0,usqr

                    if ngauss gt 1 then begin

                        muhat = ngroup * mean(xi[group]) / tausqr[k,i] + mu0[i] / usqr[i]
                        
                        muvar = 1d / (1d / usqr[i] + ngroup / tausqr[k,i])

                    endif else begin

                        muhat = ngroup * mean(xi[group]) / tausqr[k,i]

                        muvar = tausqr[k,i] / ngroup

                    endelse

                    muhat = muvar * muhat
                    
                    mu[k,i] = muhat + sqrt(muvar) * randomn(seed)
                    
                                ;get tausqr|Xi,G,mu,wsqr,nut

                    if ngauss gt 1 then begin
                        
                        nuk = ngroup + nut
                        tsqr = (nut * wsqr[i] + total( (xi[group] - mu[k,i])^2 )) / nuk
                        
                    endif else begin

                        nuk = ngroup
                        tsqr = total( (xi[group] - mu[k,i])^2 ) / nuk

                    endelse
                
                    tausqr[k,i] = tsqr * nuk / randomchi(seed, nuk)

                endif else begin

                    mu[k,i] = mu0[i] + sqrt(usqr[i]) * randomn(seed)
                    tausqr[k,i] = wsqr[i] * nut / randomchi(seed, nut)

                endelse

            endfor
                                ;get pi|G
            if ngauss eq 1 then pi[*,i] = 1d else $
              pi[*,i] = randomdir(seed, nk + 1)

        endif else begin
                                ;do metropolis-hastings updating using
                                ;approximate Gibbs sampler

            for k = 0, ngauss - 1 do begin
                
                group = where(Glabel[*,i] eq k, ngroup)
                nk[k] = ngroup

                if ngroup gt 0 then begin
                                ;get proposal for mu[k], do
                                ;approximate Gibbs sampler
                    muprop = mu[*,i]

                    muvarx = (tausqr[k,i] + mean(xvar[group]))

                    muvar = ngauss gt 1 ? 1d / (1d / usqr[i] + ngroup / muvarx) : $
                      muvarx / ngroup

                    muhat = muprop
                    
                    chisqr = randomchi(seed, 4)
                                ;draw proposal for mu from Student's t
                                ;with 4 degrees of freedom
                    muprop[k] = muhat + sqrt(muvar * 4 / chisqr) * randomn(seed)

                endif else begin
                    
                    muprop = mu[*,i]
                    muprop[k] = mu[k,i] + sqrt(usqr[i]) * randomn(seed)

                endelse

                theta = [alpha[i], beta[i], sigsqr[i]]
                
                loglik_new = loglik_mixerr( x, ygibbs, xvar, yvar, xycov, delta, theta, $
                                            pi[*,i], muprop, tausqr[*,i], Glabel[*,i] )
                logprior_new = logprior_mixerr(muprop, mu0[i], tausqr[*,i], usqr[i], wsqr[i])
                
                logpost_new = loglik_new + logprior_new
                
                accept = linmix_metro_update( logpost_new, logpost[i], seed )
                
                if accept then begin
                    
                    naccept[2+k] = naccept[2+k] + 1L
                    mu[k,i] = muprop[k]
                    logpost[i] = logpost_new
                    
                endif

                                ;get proposal for tausqr[k], do
                                ;approximate Gibbs sampler
                tsqrprop = tausqr[*,i]
 
                dof = ngroup > 1

                tsqrprop[k] = tausqr[k,i] * dof / randomchi(seed, dof)

                log_jrat = (dof + 1d) * alog(tsqrprop[k] / tausqr[k,i]) + $
                  dof / 2d * (tausqr[k,i] / tsqrprop[k] - tsqrprop[k] / tausqr[k,i])

                loglik_new = loglik_mixerr( x, ygibbs, xvar, yvar, xycov, delta, theta, $
                                            pi[*,i], mu[*,i], tsqrprop, Glabel[*,i] )
                logprior_new = logprior_mixerr(mu[*,i], mu0[i], tsqrprop, usqr[i], wsqr[i])
                
                logpost_new = loglik_new + logprior_new
                
                accept = linmix_metro_update( logpost_new, logpost[i], seed, log_jrat)
                
                if accept then begin
                    
                    naccept[2 + k + ngauss] = naccept[2 + k + ngauss] + 1L
                    tausqr[k,i] = tsqrprop[k]