<jupyter_start><jupyter_text>NGC7793 Bayesian Analysis In this notebook, we define two different models that could fit our NGC7793 data via PyMC3 and SciPy. Then we compare both of them using LOO cross validation.The data we'll use here is the temprature from the inner disk of the source (Th) and its luminosity (Luminosity_h). The same process has to be followed if we want to fit also the outer accretion disk.<jupyter_code>from IPython.display import Image

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as pm
import scipy
import scipy.stats as stats
import scipy.optimize as opt
import statsmodels.api as sm
%matplotlib inline
plt.style.use('seaborn-white')

plt.rcParams['font.family'] = 'serif'
#plt.rcParams['font.size'] = 13
#plt.rcParams['axes.labelsize'] = 10
#plt.rcParams['axes.labelweight'] = 'bold'
#plt.rcParams['axes.titlesize'] = 10
#plt.rcParams['legend.fontsize'] = 10
plt.rcParams['font.size'] = 17.0
plt.rcParams['axes.linewidth'] = 2
data=pd.read_csv('FluxData2_both.csv')
data<jupyter_output><empty_output><jupyter_text>- Modelo 1 $\propto T^4$<jupyter_code>from scipy.optimize import curve_fit
d=3.38*3.0857*10**24
data['Flux_h(cm)']=10**data['Flux_h']
data['Luminosity_h']=4*np.pi*d**2*data['Flux_h(cm)']

data['Flux_l(cm)']=10**data['Flux_l']
data['Luminosity_l']=4*np.pi*d**2*data['Flux_l(cm)']
data['Luminosity_l error(+)']=data['Flux_h']*4*np.pi*d**2 *10**(data['Flux_h']-1)* data['Flux_h error(+) ']
data['Luminosity_l error(-)']=data['Flux_h']*4*np.pi*d**2 * 10**(data['Flux_h']-1)*data['Flux_h error(-) ']

data['Luminosity_h error(+)']=data['Flux_l']*4*np.pi*d**2 * 10**(data['Flux_l']-1)*data['Flux_l error(+) ']
data['Luminosity_h error(-)']=data['Flux_l']*4*np.pi*d**2 * 10**(data['Flux_l']-1)*data['Flux_l error(-) ']

wX3_both=data.drop(11)
#Fit 1
x1=np.linspace(0.8,2.25,10000)

def func(x, A):
    return A*4*np.pi*d**2*x**4

params, covs = curve_fit(func,wX3_both['Th'],wX3_both['Luminosity_h'])
A = params[0] #stefan-boltzmann constant?
yfitx3 = A*4*np.pi*d**2*x1**4

print('σ =', A)
wX3_both['logLum_h']=np.log10(wX3_both['Luminosity_h'])
wX3_both['logTh']=np.log10(wX3_both['Th'])
logd=np.log10(d)
#Data
T=wX3_both['Th']
lum=wX3_both['logLum_h']
dis=logd
plt.scatter(np.log10(T),lum)
plt.plot(np.log10(x1),np.log10(yfitx3))
plt.title('SciPy Regression')
plt.grid()
#Model 1

with pm.Model() as model1:
    
    #Priors
    
    A1=pm.Normal(name='A1',mu=0,sigma=50)  #Constant A
    
    sigma = pm.HalfNormal("sigma", sigma=10)  #Standard deviation
        
    #Posterior
    
    Lum=np.log10(4*np.pi)+A1+2*np.log10(d)+4*np.log10(T) #Decimal logarithm luminosity formula 
    
    #Likelihood
    
    Y_obs=pm.Normal('Y_obs',Lum,sigma,observed=lum)  #lum is our luminosity data but in logarithmic scale
    
    # Prior sampling, trace definition and posterior sampling
    trace= pm.sample(draws=200000,chains=4,tune=100000,return_inferencedata=False)
import arviz as az
az.plot_trace(trace)
from arviz import plot_posterior

plot_posterior(trace);<jupyter_output>/Users/alba/my_env/lib/python3.8/site-packages/arviz/data/io_pymc3.py:96: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.
  warnings.warn(<jupyter_text>>  Regresion with bayesian model<jupyter_code>plt.figure(figsize=(7,5))
x_p=np.linspace(0.8,2.25,500)

y_p=10**(np.log10(4*np.pi)+trace['A1'].mean()+2*np.log10(d)+4*np.log10(x_p))

Lum=np.log10(4*np.pi)+A1+2*np.log10(d)+4*np.log10(T)


plt.plot(np.log10(x1),np.log10(yfitx3),label='SciPy',linewidth=4) #scipy
plt.plot(np.log10(x_p),np.log10(y_p),linestyle='dashed',label='PyMC3',linewidth=4) #bayes
plt.scatter(np.log10(T),lum,color='k',zorder=3,linewidth=2)
#plt.fill_between(np.log10(x1),np.log10(yfitx3-covs[0]),np.log10(yfitx3+covs[0]),alpha=0.2)
plt.legend(loc='upper left')
plt.xlabel('$\log_{10}(T)$')
plt.ylabel('$\log_{10}(L)$')
plt.title('$\propto T^4$ model')
plt.grid()



ax = az.plot_pair(
    trace,
    var_names=["A1", "sigma"],
    kind=["scatter", "kde"],
    kde_kwargs={"fill_last": False},
    marginals=True,
    point_estimate="median",
    figsize=(5, 5),
)

plt.show()
pm.model_to_graphviz(model1)

az.style.use("arviz-darkgrid")

az.plot_pair(
   trace, var_names=["A1", "sigma"], divergences=True, textsize=22
)
pm.model_to_graphviz(model1)<jupyter_output><empty_output><jupyter_text>-  Modelo 2  $\propto T^2$<jupyter_code>#Model 2

with pm.Model() as model2:
    
    #Priors
    
    B1=pm.Normal(name='B1',mu=0,sigma=50) #cambia mucho con el valor inicial de mu y sd
    
    sigma = pm.HalfNormal("sigma", sigma=10)
        
    #Posterior
    
    Lum=np.log10(4*np.pi)+B1+2*np.log10(d)+2*np.log10(T)
    
    #Verosimilitud
    
    Y_obs=pm.Normal('Y_obs',Lum,sigma,observed=lum)
    
    # Prior sampling, trace definition and posterior sampling
    trace1= pm.sample(draws=200000,chains=4,tune=100000,return_inferencedata=False)
from arviz import plot_posterior

plot_posterior(trace1);

#Fit 2
x2=np.linspace(0.8,2.5,10000)

def func(x, B):
    return B*4*np.pi*d**2*x**2

params, covs = curve_fit(func,wX3_both['Th'],wX3_both['Luminosity_h'])
B = params[0] #boltzmann constant
yfitx2 = B*4*np.pi*d**2*x1**2

print('B =', B)
plt.figure(figsize=(7,5))

x_p2=np.linspace(0.8,2.5,500)

y_p2=10**(np.log10(4*np.pi)+trace1['B1'].mean()+2*np.log10(d)+2*np.log10(x_p))

plt.plot(np.log10(x2),np.log10(yfitx2),label='SciPy',linewidth=4) #scipy
plt.plot(np.log10(x_p2),np.log10(y_p2),linestyle='dashed',label='PyMC3',linewidth=4) #bayes
plt.scatter(np.log10(T),lum,color='k',zorder=3,linewidth=2)
#plt.fill_between(np.log10(x1),np.log10(yfitx3-covs[0]),np.log10(yfitx3+covs[0]),alpha=0.2)
plt.legend(loc='upper left')
plt.xlabel('$\log_{10}(T)$')
plt.ylabel('$\log_{10}(L)$')
plt.title('$\propto T^2$ model')
plt.grid()



ax = az.plot_pair(
    trace1,
    var_names=["B1", "sigma"],
    kind=["scatter", "kde"],
    kde_kwargs={"fill_last": False},
    marginals=True,
    point_estimate="median",
    figsize=(5, 5),
)

plt.show()

az.plot_autocorr(trace, var_names=("A1", "sigma"))<jupyter_output>/Users/alba/my_env/lib/python3.8/site-packages/arviz/data/io_pymc3.py:96: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.
  warnings.warn(<jupyter_text>-  Comparación de los dos modelos<jupyter_code>model1.name = 'model1'
model2.name = 'model2'
df_comp_loo=pm.compare({model1.name: trace, model2.name: trace1})
df_comp_loo
plt.figure(figsize=(40,10))
az.plot_compare(df_comp_loo, insample_dev=False);
plt.grid()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18030876784153801, 'fraction_numerical': 0.053451791436061755, 'mean_word_length': 5.957446808510638, 'pattern_counts': {'":': 2, '<': 19, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '155519', 'license': 'MIT', 'n_tokens_mistral': 3166, 'n_tokens_neox': 2920, 'n_words': 469, 'path': 'Bayesian_version2.ipynb', 'repo_name': 'albarc3/NGC-7793-P9-ESAC', 'size': '6866'}