linmix_err.pro

basic median filter simulation

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;+
;   NAME:
;     LINMIX_ERR
;   PURPOSE:
;      Bayesian approach to linear regression with errors in both X and Y
;   EXPLANATION:
;     PERFORM LINEAR REGRESSION OF Y ON X WHEN THERE ARE MEASUREMENT
;     ERRORS IN BOTH VARIABLES. THE REGRESSION ASSUMES :
;
;                 ETA = ALPHA + BETA * XI + EPSILON
;                 X = XI + XERR
;                 Y = ETA + YERR
;
;
; HERE, (ALPHA, BETA) ARE THE REGRESSION COEFFICIENTS, EPSILON IS THE
; INTRINSIC RANDOM SCATTER ABOUT THE REGRESSION, XERR IS THE
; MEASUREMENT ERROR IN X, AND YERR IS THE MEASUREMENT ERROR IN
; Y. EPSILON IS ASSUMED TO BE NORMALLY-DISTRIBUTED WITH MEAN ZERO AND
; VARIANCE SIGSQR. XERR AND YERR ARE ASSUMED TO BE
; NORMALLY-DISTRIBUTED WITH MEANS EQUAL TO ZERO, VARIANCES XSIG^2 AND
; YSIG^2, RESPECTIVELY, AND COVARIANCE XYCOV. THE DISTRIBUTION OF XI
; IS MODELLED AS A MIXTURE OF NORMALS, WITH GROUP PROPORTIONS PI,
; MEAN MU, AND VARIANCE TAUSQR. BAYESIAN INFERENCE IS EMPLOYED, AND
; A STRUCTURE CONTAINING RANDOM DRAWS FROM THE POSTERIOR IS
; RETURNED. CONVERGENCE OF THE MCMC TO THE POSTERIOR IS MONITORED
; USING THE POTENTIAL SCALE REDUCTION FACTOR (RHAT, GELMAN ET
; AL.2004). IN GENERAL, WHEN RHAT < 1.1 THEN APPROXIMATE CONVERGENCE
; IS REACHED.
;
; SIMPLE NON-DETECTIONS ON Y MAY ALSO BE INCLUDED
;
; CALLING SEQUENCE:
;
;     LINMIX_ERR, X, Y, POST, XSIG=, YSIG=, XYCOV=, DELTA=, NGAUSS=, /SILENT, 
;                /METRO, MINITER= , MAXITER= 
;
;
; INPUTS :
;
;   X - THE OBSERVED INDEPENDENT VARIABLE. THIS SHOULD BE AN
;       NX-ELEMENT VECTOR.
;   Y - THE OBSERVED DEPENDENT VARIABLE. THIS SHOULD BE AN NX-ELEMENT
;       VECTOR.
;
; OPTIONAL INPUTS :
;
;   XSIG - THE 1-SIGMA MEASUREMENT ERRORS IN X, AN NX-ELEMENT VECTOR.
;   YSIG - THE 1-SIGMA MEASUREMENT ERRORS IN Y, AN NX-ELEMENT VECTOR.
;   XYCOV - THE COVARIANCE BETWEEN THE MEASUREMENT ERRORS IN X AND Y,
;           AND NX-ELEMENT VECTOR.
;   DELTA - AN NX-ELEMENT VECTOR INDICATING WHETHER A DATA POINT IS
;           CENSORED OR NOT. IF DELTA[i] = 1, THEN THE SOURCE IS
;           DETECTED, ELSE IF DELTA[i] = 0 THE SOURCE IS NOT DETECTED
;           AND Y[i] SHOULD BE AN UPPER LIMIT ON Y[i]. NOTE THAT IF
;           THERE ARE CENSORED DATA POINTS, THEN THE
;           MAXIMUM-LIKELIHOOD ESTIMATE (THETA) IS NOT VALID. THE
;           DEFAULT IS TO ASSUME ALL DATA POINTS ARE DETECTED, IE,
;           DELTA = REPLICATE(1, NX).
;   METRO - IF METRO = 1, THEN THE MARKOV CHAINS WILL BE CREATED USING
;           THE METROPOLIS-HASTINGS ALGORITHM INSTEAD OF THE GIBBS
;           SAMPLER. THIS CAN HELP THE CHAINS CONVERGE WHEN THE SAMPLE
;           SIZE IS SMALL OR IF THE MEASUREMENT ERRORS DOMINATE THE
;           SCATTER IN X AND Y.
;   SILENT - SUPPRESS TEXT OUTPUT.
;   MINITER - MINIMUM NUMBER OF ITERATIONS PERFORMED BY THE GIBBS
;             SAMPLER OR METROPOLIS-HASTINGS ALGORITHM. IN GENERAL,
;             MINITER = 5000 SHOULD BE SUFFICIENT FOR CONVERGENCE. THE
;             DEFAULT IS MINITER = 5000. THE MCMC IS STOPPED AFTER 
;             RHAT < 1.1 FOR ALL PARAMETERS OF INTEREST, AND THE
;             NUMBER OF ITERATIONS PERFORMED IS GREATER THAN MINITER.
;   MAXITER - THE MAXIMUM NUMBER OF ITERATIONS PERFORMED BY THE
;             MCMC. THE DEFAULT IS 1D5. THE MCMC IS STOPPED
;             AUTOMATICALLY AFTER MAXITER ITERATIONS.
;   NGAUSS - THE NUMBER OF GAUSSIANS TO USE IN THE MIXTURE
;            MODELLING. THE DEFAULT IS 3. IF NGAUSS = 1, THEN THE
;            PRIOR ON (MU, TAUSQR) IS ASSUMED TO BE UNIFORM.
;
; OUTPUT :
;
;    POST - A STRUCTURE CONTAINING THE RESULTS FROM THE MCMC. EACH
;           ELEMENT OF POST IS A DRAW FROM THE POSTERIOR DISTRIBUTION
;           FOR EACH OF THE PARAMETERS.
;
;             ALPHA - THE CONSTANT IN THE REGRESSION.
;             BETA - THE SLOPE OF THE REGRESSION.
;             SIGSQR - THE VARIANCE OF THE INTRINSIC SCATTER.
;             PI - THE GAUSSIAN WEIGHTS FOR THE MIXTURE MODEL.
;             MU - THE GAUSSIAN MEANS FOR THE MIXTURE MODEL.
;             TAUSQR - THE GAUSSIAN VARIANCES FOR THE MIXTURE MODEL.
;             MU0 - THE HYPERPARAMETER GIVING THE MEAN VALUE OF THE
;                   GAUSSIAN PRIOR ON MU. ONLY INCLUDED IF NGAUSS >
;                   1.
;             USQR - THE HYPERPARAMETER DESCRIBING FOR THE PRIOR
;                    VARIANCE OF THE INDIVIDUAL GAUSSIAN CENTROIDS
;                    ABOUT MU0. ONLY INCLUDED IF NGAUSS > 1.
;             WSQR - THE HYPERPARAMETER DESCRIBING THE `TYPICAL' SCALE
;                    FOR THE PRIOR ON (TAUSQR,USQR). ONLY INCLUDED IF
;                    NGAUSS > 1.
;             XIMEAN - THE MEAN OF THE DISTRIBUTION FOR THE
;                      INDEPENDENT VARIABLE, XI.
;             XISIG - THE STANDARD DEVIATION OF THE DISTRIBUTION FOR
;                     THE INDEPENDENT VARIABLE, XI.
;             CORR - THE LINEAR CORRELATION COEFFICIENT BETWEEN THE
;                    DEPENDENT AND INDEPENDENT VARIABLES, XI AND ETA.
;
; CALLED ROUTINES :
;
;    RANDOMCHI, MRANDOMN, RANDOMGAM, RANDOMDIR, MULTINOM
;
; REFERENCES :
;
;   Carroll, R.J., Roeder, K., & Wasserman, L., 1999, Flexible
;     Parametric Measurement Error Models, Biometrics, 55, 44
;
;   Kelly, B.C., 2007, Some Aspects of Measurement Error in
;     Linear Regression of Astronomical Data, ApJ, In press
;     (astro-ph/0705.2774)
;
;   Gelman, A., Carlin, J.B., Stern, H.S., & Rubin, D.B., 2004,
;     Bayesian Data Analysis, Chapman & Hall/CRC
;
; REVISION HISTORY
;
;     AUTHOR : BRANDON C. KELLY, STEWARD OBS., JULY 2006
;   - MODIFIED PRIOR ON MU0 TO BE UNIFORM OVER [MIN(X),MAX(X)] AND
;     PRIOR ON USQR TO BE UNIFORM OVER [0, 1.5 * VARIANCE(X)]. THIS
;     TENDS TO GIVE BETTER RESULTS WITH FEWER GAUSSIANS. (B.KELLY, MAY
;     2007)
; -
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;routine to compute the hyperbolic arctangent
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

function linmix_atanh, x

z = 0.5d * ( alog(1 + x) - alog(1 - x) )

return, z
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;routine to compute a robust estimate for the standard deviation of a
;data set, based on the inter-quartile range
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

function linmix_robsig, x

nx = n_elements(x)
                                ;get inter-quartile range of x
sorted = sort(x)
iqr = x[sorted[3 * nx / 4]] - x[sorted[nx / 4]]
sdev = stddev(x, /nan)
sigma = min( [sdev, iqr / 1.34] ) ;use robust estimate for sigma
if sigma eq 0 then sigma = sdev

return, sigma
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;routine to compute the log-likelihood of the data
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

function loglik_mixerr, x, y, xvar, yvar, xycov, delta, theta, pi, mu, tausqr, Glabel

alpha = theta[0]
beta = theta[1]
sigsqr = theta[2]

nx = n_elements(x)
ngauss = n_elements(pi)

Sigma11 = dblarr(nx, ngauss)
Sigma12 = dblarr(nx, ngauss)
Sigma22 = dblarr(nx, ngauss)
determ = dblarr(nx, ngauss)

for k = 0, ngauss - 1 do begin

    Sigma11[0,k] = beta^2 * tausqr[k] + sigsqr + yvar
    Sigma12[0,k] = beta * tausqr[k] + xycov
    Sigma22[0,k] = tausqr[k] + xvar

    determ[0, k] = Sigma11[*,k] * Sigma22[*,k] - Sigma12[*,k]^2

endfor

det = where(delta eq 1, ndet, comp=cens, ncomp=ncens) ;any non-detections?

loglik = dblarr(nx)

if ndet gt 0 then begin
                                ;compute contribution to
                                ;log-likelihood from the detected
                                ;sources
    for k = 0, ngauss - 1 do begin
        
        gk = where(Glabel[det] eq k, nk)

        if nk gt 0 then begin

            zsqr = (y[det[gk]] - alpha - beta * mu[k])^2 / Sigma11[det[gk],k] + $
              (x[det[gk]] - mu[k])^2 / Sigma22[det[gk],k] - $
              2d * Sigma12[det[gk],k] * (y[det[gk]] - alpha - beta * mu[k]) * $
              (x[det[gk]] - mu[k]) / (Sigma11[det[gk],k] * Sigma22[det[gk],k])
            
            corrz = Sigma12[det[gk],k] / sqrt( Sigma11[det[gk],k] * Sigma22[det[gk],k] )
            
            loglik[det[gk]] = -0.5d * alog(determ[det[gk],k]) - 0.5 * zsqr / (1d - corrz^2)

        endif

    endfor

endif

if ncens gt 0 then begin
                                ;compute contribution to the
                                ;log-likelihood from the
                                ;non-detections
    for k = 0, ngauss - 1 do begin

        gk = where(Glabel[cens] eq k, nk)

        if nk gt 0 then begin
            
            loglikx = -0.5 * alog(Sigma22[cens[gk],k]) - $
              0.5 * (x[cens[gk]] - mu[k])^2 / Sigma22[cens[gk],k]
            
                                ;conditional mean of y, given x and
                                ;G=k
            cmeany = alpha + beta * mu[k] + Sigma12[cens[gk],k] / Sigma22[cens[gk],k] * $
              (x[cens[gk]] - mu[k])
                                ;conditional variance of y, given x
                                ;and G=k
            cvary = Sigma11[cens[gk],k] - Sigma12[cens[gk],k]^2 / Sigma22[cens[gk],k]

                                ;make sure logliky is finite
            logliky = alog(gauss_pdf( (y[cens[gk]] - cmeany) / sqrt(cvary) )) > (-1d300) 
            
            loglik[cens[gk]] = loglikx + logliky
        
        endif

    endfor
    
endif

loglik = total(loglik)

return, loglik
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;routine to compute the log-prior of the data
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

function logprior_mixerr, mu, mu0, tausqr, usqr, wsqr

ngauss = n_elements(mu)

if ngauss gt 1 then begin

    logprior_mu = -0.5 * alog(usqr) - 0.5 * (mu - mu0)^2 / usqr
    logprior_mu = total(logprior_mu)
    
    logprior_tausqr = 0.5 * alog(wsqr) - 1.5 * alog(tausqr) - 0.5 * wsqr / tausqr
    logprior_tausqr = total(logprior_tausqr)

    logprior = logprior_mu + logprior_tausqr

endif else logprior = 0d ;if ngauss = 1 then uniform prior

return, logprior
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;routine to perform the Metropolis update for the scale parameter in
;the Gibbs sampler
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

function linmix_metro_update, logpost_new, logpost_old, seed, log_jrat

lograt = logpost_new - logpost_old

if n_elements(log_jrat) gt 0 then lograt = lograt + log_jrat

accept = 0

if lograt gt 0 then accept = 1 else begin
    
    u = randomu(seed)

    if alog(u) le lograt then accept = 1
    
endelse

return, accept
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;routine to acceptance rates for metropolis-hastings algorithm
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

pro linmix_metro_results, arate, ngauss

print, ''
print, 'Metropolis-Hastings Acceptance Rates:'

print, '(ALPHA, BETA) : ' + strtrim(arate[0], 1)
print, 'SIGMA^2 : ' + strtrim(arate[1], 1)
print, ''
for k = 0, ngauss - 1 do begin

    print, 'GAUSSIAN ' + strtrim(k+1,1)
    print, '   MEAN : ' + strtrim(arate[2+k], 1)
    print, '   VARIANCE : ' + strtrim(arate[2+k+ngauss], 1)

endfor

if ngauss gt 1 then begin

    print, ''
    print, 'Mu0 : ' + strtrim(arate[2+2*ngauss], 1)