<jupyter_start><jupyter_text>Read in data and define run<jupyter_code>data=pickle.load(open('/Users/cosalyk/Dropbox/visir/nirspec_fits/fztau_fluxes.p','rb'))

myrun=sf_run(data)<jupyter_output><empty_output><jupyter_text>Define prior, likelihood, posterior.Note: Posterior is prior x likelihood.  Therefore, log posterior is log prior + log likelihood.<jupyter_code>def logposterior(theta, data, sigma, myrun):

    """                                                                                       
    The natural logarithm of the joint posterior.                                             
                                                                                              
    Args:                                                                                     
        theta (list): a sample containing individual parameter values                        
        data (list): the set of data/observations                                             
        sigma (float): the standard deviation of the data points                              
        x (list): the abscissa values at which the data/model is defined                      
    """

    lp = logprior(theta) # get the prior                                                      

    # if the prior is not finite return a probability of zero (log probability of -inf)       
    if not np.isfinite(lp):
        return -np.inf

    # return the likeihood times the prior (log likelihood plus the log prior)   
    return lp + loglikelihood(theta, data, sigma, myrun)
#-------------------------    
def logprior(theta):
    """                                                                                       
    The natural logarithm of the prior probability.                                           
                                                                                              
    Args:                                                                                     
        theta (list): a sample containing individual parameter values                                                         
    """
    lp = 0.  #initialize log prior

    # unpack the model parameters from the list                                             
    logntot, temp, logomega = theta
#-------------------------    
#First parameter: logntot
    # uniform prior on log of ntot                                                                      
    lognmin = 16.  # lower range of prior                                                        
    lognmax = 24.  # upper range of prior                                                        

    # set prior to 1 (log prior to 0) if in the range and zero (log prior to -inf)
    # outside the range       
    lp = 0. if lognmin < logntot < lognmax else -np.inf
#-------------------------    
#Second parameter: temperature
    # Gaussian prior on temp, except if T<=0                                                                   
    tmu = 750.     # mean of the Gaussian prior                                                 
    tsigma = 100. # standard deviation of the Gaussian prior                                   
    lp -= 0.5*((temp - tmu)/tsigma)**2   #Add log prior due to temperature to lp due to logn.

    # set prior to zero (log prior to -inf) if T<=0       
    if (temp <= 0):
        lp = -np.inf    
#-------------------------    
#Third parameter: Omega
    #Uniform prior on Omega, but with some cutoffs
    logomegamin=-21   #0.001 AU circle at 140 pc
    logomegamax=-12   #20 AU circle at 140 pc
    lpo = 0. if logomegamin < logomega < logomegamax else -np.inf

    lp += lpo #Add log prior due to omega to lp due to temperature,logn

    return lp
#-------------------------    
def loglikelihood(theta, data, sigma, myrun):

    """                                                                                       
    The natural logarithm of the joint likelihood.                                            
                                                                                              
    Args:                                                                                     
        theta (list): a sample containing individual parameter values                        
        data (list): the set of data/observations                                             
        sigma (float): the standard deviation of the data points                              
        x (list): the abscissa values at which the data/model is defined                                       
    """

    # unpack the model parameters from the tuple                                              
    logntot, temp, logomega = theta
    omega=10**logomega
    # evaluate the model (assumes that the model is defined above)           
    md = compute_fluxes(myrun, logntot, temp, omega)
    
    # return the log likelihood  
    return -0.5*np.sum(((md - data)/sigma)**2)<jupyter_output><empty_output><jupyter_text>Run emcee<jupyter_code>Nens = 500   # number of ensemble points                                                         

#Column density: uniform prior
lognmin = 20. # lower range of prior                                                        
lognmax = 24.  # upper range of prior                                                        
lognini = np.random.uniform(lognmin, lognmax, Nens) # initial logn points                                    

#Temperature: Gaussian prior
tmu = 750.     # mean of the Gaussian prior                                                 
tsigma = 100. # standard deviation of the Gaussian prior                                   
tini = np.random.normal(tmu, tsigma, Nens) # initial t points                                    
                                      
#Omega: uniform prior    
logomegamin = -21 # lower range of prior                                                        
logomegamax = -12  # upper range of prior                                                        
logomegaini = np.random.uniform(logomegamin, logomegamax, Nens) # initial omegaini points                                    

inisamples = np.array([lognini, tini, logomegaini]).T # initial samples                                          
#100 x 2  (Nens x Nparameters)

ndims = inisamples.shape[1] # number of parameters/dimensions                                    

Nburnin = 50   # number of burn-in samples                                                      
Nsamples = 100  # number of final posterior samples                                              

# set additional args for the posterior (the data, the noise std. dev., and the abscissa)        
argslist = (myrun.lineflux, myrun.lineflux_err, myrun)

# set up the sampler                                                                             
sampler = emcee.EnsembleSampler(Nens, ndims, logposterior, args=argslist)

# pass the initial samples and total number of samples required                                  
start_time=time.time()
sampler.run_mcmc(inisamples, Nsamples+Nburnin);  
#initial start points of each walker, total number of samples
end_time=time.time()
print(Nens*(Nsamples+Nburnin), end_time-start_time)

#75000 281 seconds - Vassar machine
#20000 36 seconds - Vassar machine
#20000 73 seconds - personal machine
#4000 12.1 seconds
#2000 5.8 seconds
#200 0.5 seconds - after removing double loop<jupyter_output>75000 280.9986741542816<jupyter_text>Corner plot<jupyter_code>postsamples = sampler.chain[:, Nburnin:, :].reshape((-1, ndims))

print('Number of posterior samples is {}'.format(postsamples.shape[0]))

fig = corner.corner(postsamples, labels=[ r"$\log(n_\mathrm{tot})$",r"Temperature [K]", "Omega"])<jupyter_output>Number of posterior samples is 50000<jupyter_text>Best fit parameters with error bars<jupyter_code>from IPython.display import display, Math
labels=[ r"log(Ntot)",r"Temperature", "Omega"]
for i in range(3):
    mcmc = np.percentile(postsamples[:, i], [16, 50, 84])
    q = np.diff(mcmc)
    txt = "\mathrm{{{3}}} = {0:.3f}_{{-{1:.3f}}}^{{{2:.3f}}}"
    txt = txt.format(mcmc[1], q[0], q[1], labels[i])
    display(Math(txt))<jupyter_output><empty_output><jupyter_text>Plot of best fit<jupyter_code>bestfit=[np.percentile(postsamples[:,i],50) for i in range(3)]
modelflux=compute_fluxes(myrun, bestfit[0],bestfit[1],10**(bestfit[2]))
rotdiag=make_rotation_diagram(data,modelfluxes=modelflux)

fig=plt.figure()
ax1=fig.add_subplot(111)
ax1.plot(rotdiag['x'],rotdiag['y'],'ro',label='Data')
ax1.plot(rotdiag['x'],rotdiag['modely'],'bo',label='Model')
ax1.plot(rotdiag['x'],rotdiag['modely'],'b',label='_None_')
ax1.set_xlabel('Upper Level Energy [Kelvin]')
ax1.legend()<jupyter_output><empty_output><jupyter_text>Plot of some random samples<jupyter_code>fig=plt.figure()
ax1=fig.add_subplot(111)
inds = np.random.randint(len(postsamples), size=100)
for ind in inds:
    sample = postsamples[ind]
    modelflux=compute_fluxes(myrun, sample[0],sample[1],10.**sample[2])
    rotdiag=make_rotation_diagram(data,modelfluxes=modelflux)
    ax1.plot(rotdiag['x'],rotdiag['modely'],alpha=0.2)
ax1.set_xlabel('Upper Level Energy [Kelvin]')

bestfit=[np.percentile(postsamples[:,i],50) for i in range(3)]
modelflux=compute_fluxes(myrun, bestfit[0],bestfit[1],10**(bestfit[2]))
rotdiag=make_rotation_diagram(data,modelfluxes=modelflux)
ax1.plot(rotdiag['x'],rotdiag['y'],'ro',label='Data')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09808890296695175, 'fraction_numerical': 0.020694752402069475, 'mean_word_length': 1.3816947447824994, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 85, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '567783', 'license': 'MIT-0', 'n_tokens_mistral': 2548, 'n_tokens_neox': 2328, 'n_words': 674, 'path': '.ipynb_checkpoints/slab_fitter-checkpoint.ipynb', 'repo_name': 'csalyk/slab_fitter', 'size': '9471'}