linmix_err.pro

basic median filter simulation

                    logpost[i] = logpost_new
                    
                endif
                
            endfor
                                ;get pi|G, can do exact Gibbs sampler
                                ;for this
            if ngauss eq 1 then pi[*,i] = 1d else $
              pi[*,i] = randomdir(seed, nk + 1)

        endelse
        
;finally, update parameters for prior distribution, only do this if
;more than one gaussian
        
        if ngauss gt 1 then begin

            if gibbs then begin
                
                repeat mu0[i] = mean(mu[*,i]) + sqrt(usqr[i] / ngauss) * randomn(seed) $
                  until (mu0[i] gt mu0min) and (mu0[i] lt mu0max)

            endif else begin
                
                loglik = loglik_mixerr( x, ygibbs, xvar, yvar, xycov, delta, theta, $
                                        pi[*,i], mu[*,i], tausqr[*,i], Glabel[*,i] )
                
                muprop = mu0[i] + sqrt(jvar_mu0) * randomn(seed)

                if muprop gt mu0min and muprop lt mu0max then begin

                    logprior_old = logprior_mixerr(mu[*,i], mu0[i], tausqr[*,i], usqr[i], wsqr[i])
                    logprior_new = logprior_mixerr(mu[*,i], muprop, tausqr[*,i], usqr[i], wsqr[i])
                    
                    logpost_new = loglik + logprior_new
                    logpost_old = loglik + logprior_old
                    
                    accept = linmix_metro_update( logpost_new, logpost_old, seed )
                    
                    if accept then begin
                        
                        naccept[2 + 2 * ngauss] = naccept[2 + 2 * ngauss] + 1L
                        mu0[i] = muprop
                        logpost[i] = loglik + logprior_new
                        
                    endif
                
                endif

            endelse

            if gibbs then begin
                
                nu = ngauss + nuu
                usqr0 = (nuu * wsqr[i] + total( (mu[*,i] - mu0[i])^2 )) / nu
                
                repeat usqr[i] = usqr0 * nu / randomchi(seed, nu) $
                  until usqr[i] le umax
                
            endif else begin
                                ;do metropolis update

                log_usqr = alog(usqr[i]) + sqrt(jvar_usqr) * randomn(seed)
                usqr0 = exp(log_usqr)
                
                if usqr0 le umax then begin

                    logprior_old = logprior_mixerr(mu[*,i], mu0[i], tausqr[*,i], usqr[i], wsqr[i])
                    
                    logpost[i] = loglik + logprior_old ;update posterior after gibbs step for mu0
                    
                    logprior_new = logprior_mixerr(mu[*,i], mu0[i], tausqr[*,i], usqr0, wsqr[i])
                    
                    logpost_new = loglik + logprior_new
                    logpost_old = loglik + logprior_old
                    
                    log_jrat = log_usqr - alog(usqr[i])
                    
                    accept = linmix_metro_update( logpost_new, logpost_old, seed, log_jrat )
                    
                    if accept then begin
                        
                        naccept[3 + 2 * ngauss] = naccept[3 + 2 * ngauss] + 1L
                        usqr[i] = usqr0
                        logpost[i] = loglik + logprior_new
                        
                    endif

                endif
                
            endelse
            
            if gibbs then begin
                
                alphaw = ngauss * nut / 2d + 1
                betaw = 0.5 * nut * total(1d / tausqr[*,i])
                
                wsqr[i] = randomgam(seed, alphaw, betaw)
                
            endif else begin
                
                log_wsqr = alog(wsqr[i]) + sqrt(jvar_wsqr) * randomn(seed)
                wsqr0 = exp(log_wsqr)
                
                logprior_old = logprior_mixerr(mu[*,i], mu0[i], tausqr[*,i], usqr[i], wsqr[i])
                logprior_new = logprior_mixerr(mu[*,i], mu0[i], tausqr[*,i], usqr[i], wsqr0)
                
                logpost_new = loglik + logprior_new + log_wsqr
                logpost_old = loglik + logprior_old + alog(wsqr[i])
                
                accept = linmix_metro_update( logpost_new, logpost_old, seed )
                
                if accept then begin
                    
                    naccept[4 + 2 * ngauss] = naccept[4 + 2 * ngauss] + 1L
                    wsqr[i] = wsqr0
                    logpost[i] = loglik + logprior_new
                    
                endif
                
            endelse

        endif

    endfor

                                ;save Markov Chains
    if iter eq 0 then begin
        
        alphag = alpha
        betag = beta
        sigsqrg = sigsqr
        
        pig = pi
        mug = mu
        tausqrg = tausqr

        if ngauss gt 1 then begin

            mu0g = mu0
            usqrg = usqr
            wsqrg = wsqr
        
        endif

        if metro then logpostg = logpost

    endif else begin
        
        alphag = [alphag, alpha]
        betag = [betag, beta]
        sigsqrg = [sigsqrg, sigsqr]

        pig = [[pig], [pi]]
        mug = [[mug], [mu]]
        tausqrg = [[tausqrg], [tausqr]]
        
        if ngauss gt 1 then begin
        
            mu0g = [mu0g, mu0]
            usqrg = [usqrg, usqr]
            wsqrg = [wsqrg, wsqr]

        endif
        
        if metro then logpostg = [logpostg, logpost]

    endelse
        
    iter = iter + 1L
    
;check for convergence
    
    if iter ge 4 and iter eq checkiter and not burnin then begin

        if not silent and metro then linmix_metro_results, $
          float(naccept) / (nchains * iter), ngauss

        Bvar = dblarr(npar)     ;between-chain variance
        Wvar = dblarr(npar)     ;within-chain variance
        
        psi = dblarr(iter, nchains, npar)
        
        psi[*,*,0] = transpose(reform(alphag, nchains, iter))
        psi[*,*,1] = transpose(reform(betag, nchains, iter))
        psi[*,*,2] = transpose(reform(sigsqrg, nchains, iter))

        pig2 = reform(pig, ngauss, nchains, iter)
        mug2 = reform(mug, ngauss, nchains, iter)
        tausqrg2 = reform(tausqrg, ngauss, nchains, iter)

        psi[*,*,3] = transpose( total(pig2 * mug2, 1) ) ;mean of xi
                                ;variance of xi
        psi[*,*,4] = transpose( total(pig2 * (tausqrg2 + mug2^2), 1) ) - psi[*,*,3]^2
                                ;linear correlation coefficient
                                ;between xi and eta
        psi[*,*,5] = psi[*,*,1] * sqrt(psi[*,*,4] / (psi[*,*,1]^2 * psi[*,*,4] + psi[*,*,2]))
                                ;do normalizing transforms before
                                ;monitering convergence
        psi[*,*,2] = alog(psi[*,*,2])
        psi[*,*,4] = alog(psi[*,*,4])
        psi[*,*,5] = linmix_atanh(psi[*,*,5])

        psi = psi[iter/2:*,*,*] ;discard first half of MCMC

        ndraw = iter / 2
                                ;calculate between- and within-sequence
                                ;                  variances
        for j = 0, npar - 1 do begin
            
            psibarj = total( psi[*,*,j], 1 ) / ndraw
            psibar = mean(psibarj)
            
            sjsqr = 0d
            for i = 0, nchains - 1 do $
              sjsqr = sjsqr + total( (psi[*, i, j] - psibarj[i])^2 ) / (ndraw - 1.0)
            
            Bvar[j] = ndraw / (nchains - 1.0) * total( (psibarj - psibar)^2 )
            Wvar[j] = sjsqr / nchains
            
        endfor
        
        varplus = (1.0 - 1d / ndraw) * Wvar + Bvar / ndraw
        Rhat = sqrt( varplus / Wvar ) ;potential variance scale reduction factor
        
        if total( (Rhat le 1.1) ) eq npar and iter ge miniter then convergence = 1 $
        else if iter ge maxiter then convergence = 1 else begin
            
            if not silent then begin
                print, 'Iteration: ', iter
                print, 'Rhat Values for ALPHA, BETA, log(SIGMA^2), mean(XI), ' + $
                  'log(variance(XI), atanh(corr(XI,ETA)) ): '
                print, Rhat
            endif
            
            checkiter = checkiter + 100L
            
        endelse

    endif
    
    if (burnin) and (iter eq burniter) then begin
;still doing burn-in stage, get new estimates for jumping kernel
;parameters
        
        jvar_ssqr = linmix_robsig( alog(sigsqrg) )^2

                                ;now modify covariance matrix for
                                ;coefficient jumping kernel
        coefg = [[alphag], [betag]]
        
        jvar_coef = correlate( transpose(coefg), /covar)

        if ngauss gt 1 then begin

            jvar_mu0 = linmix_robsig(mu0g)^2 * 2.4^2

            jvar_usqr = linmix_robsig( alog(usqrg) )^2 * 2.4^2
    
            jvar_wsqr = linmix_robsig( alog(wsqrg) )^2 * 2.4^2

        endif
        
        if iter eq burnstop then burnin = 0
        
        if not silent and not burnin then begin
            
            print, 'Burn-in Complete'
            iter = 0L
            
        endif
        
        naccept = lonarr(5 + 2 * ngauss)
        burniter = burniter + 250L
        
    endif
    
endrep until convergence

ndraw = iter * nchains / 2

;save posterior draws in a structure

if ngauss gt 1 then begin

    post = {alpha:0d, beta:0d, sigsqr:0d, pi:dblarr(ngauss), mu:dblarr(ngauss), $
            tausqr:dblarr(ngauss), mu0:0d, usqr:0d, wsqr:0d, ximean:0d, xisig:0d, $
            corr:0d}

endif else begin
    
    post = {alpha:0d, beta:0d, sigsqr:0d, pi:dblarr(ngauss), mu:dblarr(ngauss), $
            tausqr:dblarr(ngauss), ximean:0d, xisig:0d, corr:0d}

endelse

post = replicate(post, ndraw)

post.alpha = alphag[(iter*nchains+1)/2:*]
post.beta = betag[(iter*nchains+1)/2:*]
post.sigsqr = sigsqrg[(iter*nchains+1)/2:*]
post.pi = pig[*,(iter*nchains+1)/2:*]
post.mu = mug[*,(iter*nchains+1)/2:*]
post.tausqr = tausqrg[*,(iter*nchains+1)/2:*]

if ngauss gt 1 then begin

    post.mu0 = mu0g[(iter*nchains+1)/2:*]
    post.usqr = usqrg[(iter*nchains+1)/2:*]
    post.wsqr = wsqrg[(iter*nchains+1)/2:*]

endif

post.ximean = total(post.pi * post.mu, 1) ;mean of xi   
post.xisig = total(post.pi * (post.tausqr + post.mu^2), 1) - post.ximean^2
post.xisig = sqrt(post.xisig)   ;standard deviation of xi

                                ;get linear correlation coefficient
                                ;between xi and eta
post.corr = post.beta * post.xisig / sqrt(post.beta^2 * post.xisig^2 + post.sigsqr)

return
end