sedfitting/fitsedcombine_mc.py at master · bjweiner/sedfitting · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103

# Derived from fitsedfamily_mc.py

# Given an SED and some error bars, generate N_mc monte carlo
# realizations of the SED, fit models to each realization with the
# fitsedcombine approach - fitting a combination and trying to make it
# sparse -
# store the predicted best fit coefficients, and the chisq,
# which will typically be pretty low

# Then compute, over the realizations, the:
#  - mean and rms of each coefficient
#  - mean and rms of the sum of the coeffs, which should be total mass

# return a list that contains statistics of the mc fits

# bjw 5/31/2016

import numpy as np
import matplotlib.pyplot as plt
import fitsedcombine

def fitsedcombine_mc(sedstruct, z, wobs, fobs, efobs, filtstruct, nmonte, penalize=0.0, initguess=0, makeplot=0, logplot=0):

    nsed = len(sedstruct)
    nobs = len(wobs)
    nbestarray = []
    bestcoeffsarray = np.zeros([nmonte,nsed])
    meancoeffarray = np.zeros(nsed)
    rmscoeffarray = np.zeros(nsed)
    sumcoeffarray = np.zeros(nmonte)
    chisqbestarray = np.zeros(nmonte)
    nzerocoeffarray = np.zeros(nmonte)

    for i in range(nmonte):
        errs_in_sigma = np.random.normal(0.0,1.0,nobs)
        fluxtmp = fobs + efobs * errs_in_sigma
        fitcoeffs, fiterrors, chisq1 = fitsedcombine.fitsedcombine(sedstruct, z, wobs, fluxtmp, efobs, filtstruct, makeplot=0, penalize=penalize, initguess=initguess)
        bestcoeffsarray[i,0:] = fitcoeffs
        chisqbestarray[i] = chisq1
        sumcoeffarray[i] = np.sum(fitcoeffs)
        nzerocoeffarray[i] = sum(fitcoeffs < 1.0e-4)

    for j in range(nsed):
        meancoeffarray[j] = np.mean(bestcoeffsarray[0:,j])
        rmscoeffarray[j] = np.std(bestcoeffsarray[0:,j],ddof=1)

    meansumcoeff = np.mean(sumcoeffarray)
    rmssumcoeff = np.std(sumcoeffarray,ddof=1)
    meannzerocoeff = np.mean(nzerocoeffarray)
    rmsnzerocoeff = np.std(nzerocoeffarray,ddof=1)

    print "mean fit coeffs per SED: ", meancoeffarray
    print "rms of fit coeffs: ", rmscoeffarray
    print "mean, rms of sum of coeffs: ", meansumcoeff, rmssumcoeff
    print "mean, rms of num of zero coeffs: ", meannzerocoeff, rmsnzerocoeff

    chisqbest_mean = np.mean(chisqbestarray)
    chisqbest_rms = np.std(chisqbestarray, ddof=1)

    print "avg, rms of chisq(bestfit): ",chisqbest_mean,chisqbest_rms

    # make a plot?  eg a histogram of the values of the coeffs
    if makeplot == 1:
        plt.clf()
        # set bins?
        colorarray = ['blue','cyan','green','red','magenta','black']
        for j in range(nsed):
            # Only make plot if the SED had some non-zero coeffs
            if np.max(bestcoeffsarray[0:,j]) > 1.0e-6:
                if logplot ==1:
                    xcoeff = np.log10(bestcoeffsarray[0:,j])
                else:
                    xcoeff = bestcoeffsarray[0:,j]
                    icolor = j % len(colorarray)
                    histvals1, binvals1, patches1 = plt.hist(xcoeff, color=colorarray[icolor], histtype='step')
        if logplot==1:
            fig = plt.xlabel('log M / Msun')
        else:
            fig = plt.xlabel('M / Msun')
        fig = plt.ylabel('number of realizations')
        # plt.figlegend( (patches1[0], patches2[0]), ('point estimates', 'expect values'), 'upper right')
        plt.show()
        plt.clf()
        if logplot==1:
            xsumcoeff = np.log10(sumcoeffarray)
            fig = plt.xlabel('log M / Msun')
        else:
            xsumcoeff = sumcoeffarray
            fig = plt.xlabel('Mass M')
        fig = plt.ylabel('number of realizations')
        histvals2, binvals2, patches2 = plt.hist(xsumcoeff, color='black', histtype='step')
        plt.show()

    # first list is quantities in units of the normalization, eg Msun
    result_coeffs = {'meansum':meansumcoeff, 'rmssum':rmssumcoeff, 'meancoeffs':meancoeffarray, 'rmscoeffs':rmscoeffarray, 'meannzero':meannzerocoeff, 'rmsnzero':rmsnzerocoeff }
    # , 'nzerocoeffs':nzerocoeffarray }
    # second list has to do with fit quality
    result_probs = [chisqbest_mean, chisqbest_rms]
    # third dict is values for each realization
    mc_fits = {'bestcoeffs': bestcoeffsarray, 'sumcoeff': sumcoeffarray}

    return result_coeffs, result_probs, mc_fits