linmix_err.pro

basic median filter simulation

    print, 'u^2 : ' + strtrim(arate[3+2*ngauss], 1)
    print, 'w^2 : ' + strtrim(arate[4+2*ngauss], 1)

endif

print, ''

return
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;                                                                     ;
;                             MAIN ROUTINE                            ;
;                                                                     ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

pro linmix_err, x, y, post, xsig=xsig, ysig=ysig, xycov=xycov, delta=delta, $
                ngauss=ngauss, metro=metro, silent=silent, miniter=miniter, $
                maxiter=maxiter

if n_params() lt 3 then begin

    print, 'Syntax- LINMIX_ERR, X, Y, POST, XSIG=XSIG, YSIG=YSIG, XYCOV=XYCOV,'
    print, '                    DELTA=DELTA, NGAUSS=NGAUSS, /SILENT, /METRO, '
    print, '                    MINITER=MINITER, MAXITER=MAXITER'
    return

endif

;check inputs and setup defaults

nx = n_elements(x)
if n_elements(y) ne nx then begin
    print, 'Y and X must have the same size.'
    return
endif

if n_elements(xsig) eq 0 and n_elements(ysig) eq 0 then begin
    print, 'Must supply at least one of XSIG or YSIG.'
    return
endif

if n_elements(xsig) eq 0 then begin
    xsig = dblarr(nx)
    xycov = dblarr(nx)
endif
if n_elements(ysig) eq 0 then begin
    ysig = dblarr(nx)
    xycov = dblarr(nx)
endif
if n_elements(xycov) eq 0 then xycov = dblarr(nx)

if n_elements(xsig) ne nx then begin
    print, 'XSIG and X must have the same size.'
    return
endif
if n_elements(ysig) ne nx then begin
    print, 'YSIG and X must have the same size.'
    return
endif
if n_elements(xycov) ne nx then begin
    print, 'XYCOV and X must have the same size.'
    return
endif

if n_elements(delta) eq 0 then delta = replicate(1, nx)
if n_elements(delta) ne nx then begin
    print, 'DELTA and X must have the same size.'
    return
endif

bad = where(finite(x) eq 0 or finite(y) eq 0 or finite(xsig) eq 0 or $
            finite(ysig) eq 0 or finite(xycov) eq 0, nbad)

if nbad gt 0 then begin
    print, 'Non-finite input detected.'
    return
endif

det = where(delta eq 1, ndet, comp=cens, ncomp=ncens) ;get detected data points

if ncens gt 0 then begin

    cens_noerr = where(ysig[cens] eq 0, ncens_noerr)
    if ncens_noerr gt 0 then begin
        print, 'NON-DETECTIONS FOR Y MUST HAVE NON-ZERO MEASUREMENT ERROR VARIANCE.'
        return
    endif

endif

 ;find data points without measurement error
xnoerr = where(xsig eq 0, nxnoerr, comp=xerr, ncomp=nxerr)
ynoerr = where(ysig eq 0, nynoerr, comp=yerr, ncomp=nyerr)

if nxerr gt 0 then ynoerr2 = where(ysig[xerr] eq 0, nynoerr2) else nynoerr2 = 0L
if nyerr gt 0 then xnoerr2 = where(xsig[yerr] eq 0, nxnoerr2) else nxnoerr2 = 0L

xvar = xsig^2
yvar = ysig^2
xycorr = xycov / (xsig * ysig)
if nxnoerr gt 0 then xycorr[xnoerr] = 0d
if nynoerr gt 0 then xycorr[ynoerr] = 0d

if not keyword_set(metro) then metro = 0
if metro then gibbs = 0 else gibbs = 1
if not keyword_set(silent) then silent = 0
if n_elements(ngauss) eq 0 then ngauss = 3

if ngauss le 0 then begin
    print, 'NGAUSS must be at least 1.'
    return
endif

if n_elements(miniter) eq 0 then miniter = 5000L ;minimum number of iterations that the 
                                                 ;Markov Chain must perform
if n_elements(maxiter) eq 0 then maxiter = 100000L ;maximum number of iterations that the 
                                                   ;Markov Chain will perform

;; perform MCMC

nchains = 4                    ;number of markov chains
checkiter = 100               ;check for convergence every 100 iterations
iter = 0L

;use BCES estimator for initial guess of theta = (alpha, beta, sigsqr)
beta = ( correlate(x, y, /covar) - mean(xycov) ) / $
  ( variance(x) - mean(xvar) )
alpha = mean(y) - beta * mean(x)

sigsqr = variance(y) - mean(yvar) - beta * (correlate(x,y, /covar) - mean(xycov))
sigsqr = sigsqr > 0.05 * variance(y - alpha - beta * x)

                                ;get initial guess of mixture
                                ;parameters prior
mu0 = median(x)
wsqr = variance(x) - median(xvar)
wsqr = wsqr > 0.01 * variance(x)

;now get MCMC starting values dispersed around these initial guesses

Xmat = [[replicate(1d, nx)], [x]]
Vcoef = invert( Xmat ## transpose(Xmat), /double ) * sigsqr

coef = mrandomn(seed, Vcoef, nchains)
chisqr = randomchi(seed, 4, nchains)

;randomly disperse starting values for (alpha,beta) from a
;multivariate students-t distribution with 4 degrees of freedom
alphag = alpha + coef[*,0] * sqrt(4d / chisqr)
betag = beta + coef[*,1] * sqrt(4d / chisqr)

                                ;draw sigsqr from an Inverse scaled
                                ;chi-square density
sigsqrg = sigsqr * (nx / 2) / randomchi(seed, nx / 2, nchains)

;get starting values for the mixture parameters, first do prior
;parameters
   
                                ;mu0 is the global mean

mu0min = min(x) ;prior for mu0 is uniform over mu0min < mu0 < mu0max
mu0max = max(x)

repeat begin
    
    mu0g = mu0 + sqrt(variance(x) / nx) * randomn(seed, nchains) / $
      sqrt(4d / randomchi(seed, 4, nchains))

    pass = where(mu0g gt mu0min and mu0g lt mu0max, npass)

endrep until npass eq nchains

                                ;wsqr is the global scale
wsqrg = wsqr * (nx / 2) / randomchi(seed, nx / 2, nchains)

usqrg = replicate(variance(x) / 2d, nchains)

;now get starting values for mixture parameters

tausqrg = dblarr(ngauss, nchains) ;initial group variances
for k = 0, ngauss - 1 do tausqrg[k,*] = 0.5 * wsqrg * 4 / $
  randomchi(seed, 4, nchains)

mug = dblarr(ngauss, nchains)   ;initial group means
for k = 0, ngauss - 1 do mug[k,*] = mu0g + sqrt(wsqrg) * randomn(seed, nchains)

;get initial group proportions and group labels

pig = dblarr(ngauss, nchains)
Glabel = intarr(nx, nchains)

if ngauss eq 1 then Glabel = intarr(nx, nchains) else begin

    for i = 0, nchains - 1 do begin
        
        for j = 0, nx - 1 do begin
                                ;classify sources to closest centroid
            dist = abs(mug[*,i] - x[j])
            mindist = min(dist, minind)
            
            pig[minind,i] = pig[minind,i] + 1
            Glabel[j,i] = minind
            
        endfor

    endfor

endelse
                                ;get initial values for pi from a
                                ;dirichlet distribution, with
                                ;parameters based on initial class
                                ;occupancies
if ngauss eq 1 then pig = transpose(replicate(1d, nchains)) else $
  for i = 0, nchains - 1 do pig[*,i] = randomdir(seed, pig[*,i] + 1)

alpha = alphag
beta = betag
sigsqr = sigsqrg
mu = mug
tausqr = tausqrg
pi = pig
mu0 = mu0g
wsqr = wsqrg
usqr = usqrg

eta = dblarr(nx, nchains)
for i = 0, nchains - 1 do eta[*,i] = y ;initial values for eta

nut = 1 ;degrees of freedom for the prior on tausqr
nuu = 1 ;degrees of freedom for the prior on usqr

;number of parameters to moniter convergence on
npar = 6

if metro then begin
;get initial variances for the jumping kernels

    jvar_coef = Vcoef
    log_ssqr = alog( sigsqr[0] * nx / randomchi(seed, nx, 1000) )
    jvar_ssqr = variance(log_ssqr) ;get variance of the jumping density
                                   ;for sigsqr

                                ;get variances for prior variance
                                ;parameters
    jvar_mu0 = variance(x) / ngauss
    jvar_wsqr = variance( alog(variance(x) * nx / randomchi(seed, nx, 1000)) )
    jvar_usqr = jvar_wsqr

    naccept = lonarr(5 + 2 * ngauss)
    
    logpost = dblarr(nchains)
                                ;get initial values of the
                                ;log-posterior
    for i = 0, nchains - 1 do begin
        
        theta = [alpha[i], beta[i], sigsqr[i]]
        loglik = loglik_mixerr( x, y, 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
        
    endfor
    
endif

convergence = 0

;stop burn-in phase after BURNSTOP iterations if doing
;Metropolis-Hastings jumps, update jumping kernels every BURNITER
;iterations

burnin = metro ? 1 : 0
burniter = 250
burnstop = 500 < (miniter / 2 > 100)
                                ;start Markov Chains
if not silent then print, 'Simulating Markov Chains...'
if not silent and metro then print, 'Doing Burn-in First...'

ygibbs = y
xi = x
umax = 1.5 * variance(x) ;prior for usqr is uniform over 0 < usqr < umax

if metro then begin
                                ;define arrays now so we don't have to
                                ;create them every MCMC iteration
    Sigma11 = dblarr(nx, ngauss)
    Sigma12 = dblarr(nx, ngauss)
    Sigma22 = dblarr(nx, ngauss)
    determ = dblarr(nx, ngauss)

endif

gamma = dblarr(nx, ngauss)
nk = fltarr(ngauss)

repeat begin
    
    for i = 0, nchains - 1 do begin ;do markov chains one at-a-time
        
        if gibbs then begin
            
            if ncens gt 0 then begin
                                ;first get new values of censored y
                for j = 0, ncens - 1 do begin
                    
                    next = 0
                    repeat ygibbs[cens[j]] = eta[cens[j],i] + $
                      sqrt(yvar[cens[j]]) * randomn(seed) $
                      until ygibbs[cens[j]] le y[cens[j]]
                    
                endfor
                
            endif
            
;need to get new values of Xi and Eta for Gibbs sampler
            
            if nxerr gt 0 then begin
                                ;first draw Xi|theta,x,y,G,mu,tausqr
                xixy = x[xerr] + xycov[xerr] / yvar[xerr] * (eta[xerr,i] - ygibbs[xerr])
                if nynoerr2 gt 0 then xixy[ynoerr2] = x[xerr[ynoerr2]]
                xixyvar = xvar[xerr] * (1 - xycorr[xerr]^2)
                
                for k = 0, ngauss - 1 do begin ;do one gaussian at-a-time
                    
                    group = where(Glabel[xerr,i] eq k, ngroup)