; Gaussian probability FUNCTION _gauss_prob, x, p, nx, np, funa = funa prob = 0D norm = 0D ; Norm FOR igauss = 0, funa.ngauss-1 DO norm = norm+p[3*igauss+2, *] norm = transpose(rebin(reform(norm), np, nx)) FOR igauss = 0, funa.ngauss-1 DO BEGIN p0 = transpose(rebin(reform(p[3*igauss, *]), np, nx)) p1 = transpose(rebin(reform(p[3*igauss+1, *]), np, nx)) p2 = transpose(rebin(reform(p[3*igauss+2, *]), np, nx))/norm prob = prob+(p2/sqrt(2*!pi*(p1^2)))*exp(-(x-p0)^2/(2*(p1^2))) ENDFOR return, prob END ; Likelihood FUNCTION _gauss_ml, p, np, funa = funa nx = n_elements(funa.x) xx = rebin(funa.x, nx, np) prob = _gauss_prob(xx, p, nx, np, funa = funa) logl = -total(alog(prob), 1) return, logl END ; SETUP binsize = 0.3 xr = [-3, 16] nmod = 1000 xmod = xr[0]+(xr[1]-xr[0])*findgen(nmod)/nmod col0 = 180 th = 5 ; DATA ngauss = 4 x = 0 means = [2, 5, 9, 12] sigmas = [1, 3, 0.5, 1] ns = [90, 140, 170, 150] FOR igauss = 0, ngauss-1 DO x = [x, means[igauss]+sigmas[igauss]*randomn(seed, ns[igauss])] x = x[1:*] n = n_elements(x) ; SOLBER call ntry = 7 probs = fltarr(ntry, n) sols = fltarr(ntry, 3*ntry) FOR itry = 0, ntry-1 DO BEGIN ngauss = itry+1 funa = {x:x, ngauss:ngauss} lim = reform(rebin(reform([[1e-1, 15], [1e-1, 5.], [0, 1]], 6), 6, ngauss), 2, 3*ngauss) sol = solber('_gauss_ml', 3*funa.ngauss, funa = funa, lim = lim, npop = 200, plot_flag = 2, mutrate = 0.05, mrate_max = 0.15, cross = 0.6, nrep = 0.5, ngen_max = 200, term_fit = 0, status = status, ngen_tot = ngen_tot, gfit_best = gfit_best, save_gen = save_gen, save_gfit = save_gfit) print, gfit_best print, sol ; Normalize weights norm = 0 FOR igauss = 0, funa.ngauss-1 DO norm = norm+sol[3*igauss+2] FOR igauss = 0, funa.ngauss-1 DO sol[3*igauss+2] = sol[3*igauss+2]/norm ; Save solution sols[itry, 0:3*itry+2] = sol ; Probability prob = _gauss_prob(x, sol, n, ngauss, funa = funa) probs[itry, *] = prob ENDFOR ; Choose Solber solution ngauss = 4 funa = {ngauss:ngauss} sol = reform(sols[ngauss-1, 0:3*ngauss-1]) ; Original model p = fltarr(3, ngauss) p[0, *] = means p[1, *] = sigmas p[2, *] = float(ns)/n p = reform(p, 3*ngauss) goris = fltarr(funa.ngauss, nmod) FOR igauss = 0, funa.ngauss-1 DO goris[igauss, *] = p[3*igauss+2]*binsize*n*_gauss_prob(xmod, p[3*igauss:3*igauss+2], nmod, 1, funa = {ngauss:1}) gori = binsize*n*_gauss_prob(xmod, p, nmod, ngauss, funa = funa) ; Histogram h = histogram(x, binsize = binsize, locations = hbins) hbins = hbins+0.5*binsize ; Solber models gmods = fltarr(funa.ngauss, nmod) FOR igauss = 0, funa.ngauss-1 DO gmods[igauss, *] = sol[3*igauss+2]*binsize*n*_gauss_prob(xmod, sol[3*igauss:3*igauss+2], nmod, 1, funa = {ngauss:1}) ; PLOT window, 1, retain = 2, xs = 950, ys = 180 !p.multi = [0, 4, 1] !y.margin = [3, 2] !x.margin = [6, 1] !p.charsize = 1.8 ; Original yr = [0, 1.2*max(h)] plot, [0], [0], xr = xr, xs = 1, yr = yr, ys = 1, backgr = 255, color = 0, title = 'Model + data' oplot, xmod, gori, color = col0, thick = th oplot, hbins, h, color = 0, psym = 10 FOR igauss = 0, funa.ngauss-1 DO oplot, xmod, goris[igauss, *], color = 0 ; Solber model plot, [0], [0], xr = xr, xs = 1, yr = yr, ys = 1, backgr = 255, color = 0, title = 'Model + data + solber solution' oplot, xmod, gori, color = col0, thick = th oplot, hbins, h, color = 0, psym = 10 FOR igauss = 0, funa.ngauss-1 DO oplot, xmod, gmods[igauss, *], color = 0 gtot = total(gmods, 1) oplot, xmod, gtot, color = 0, thick = 2 ; CDF plot, [0], [0], xr = xr, xs = 1, yr = [0, 1], ys = 1, color = 0, title = 'CDFs' oplot, xmod, total(gori, /cumulative)/total(gori), color = col0, thick = th oplot, hbins, total(h, /cumulative)/n, color = 0 oplot, xmod, total(gtot, /cumulative)/total(gtot), color = 0, thick = 2 ; Likelihood ratio lrat = fltarr(ntry-1) FOR itry = 0, ntry-2 DO lrat[itry] = -2*total(alog(probs[itry, *]/probs[itry+1, *])) plot, indgen(ntry-1)+1, lrat, xtitle = 'N gaussians', color = 0, title = 'Likelihood ratio test' oplot, [1, ntry-1], 3*[1, 1], linestyle = 1, color = 0 xyouts, 1.2, max(lrat), '-2 ln (L!Di!N/L!Di+1!N)', color = 0, charsize = 1.2 !p.multi = 0 ; Image save snapshot = TVRD(True = 1) write_png, 'solber_gaussmixture.png', snapshot END