想问下，我这个代码哪里有问题吗？能跑但好像没办法正确计算：
```python
import emcee
import corner
import numpy as np
import afterglowpy as grb
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from multiprocessing import Pool
import cal
import os
# 读取txt数据文件并整理数据
t_x, f_x, e1_x, e2_x = np.loadtxt(‘Xray.txt', usecols=(0, 1, 2, 3), unpack=True)
t_r, f_r, e_r = np.loadtxt('radio.txt', usecols=(0, 2, 3), unpack=True)
t_o, f_o, e_o = np.loadtxt('optical.txt', usecols=(0, 1, 2), unpack=True)
t_all = np.concatenate((t_r, t_o, t_x), axis=None)
y = np.concatenate((f_r, f_o, f_x), axis=None)
yerr_h = np.concatenate((e_r / 2.0, e_o / 2.0, e1_x), axis=None)
yerr_l = np.concatenate((e_r / 2.0, e_o / 2.0, -e2_x ), axis=None)
# 每个波段所占个数
num_r = len(t_r)
num_o = len(t_o)
num_x = len(t_x)
# 固定参数：
z = 1.0
################################################
# 模型函数
def flux_Model(t, log_E_e, log_E_b, log_dNe, log_E_iso_w, log_E_iso_n, log_Eta_0_w, log_Eta_0_n, theta_j_w, theta_j_n, theta_v_fj, theta_v_rj):
global z
t1 = np.sort(t)
# 计算flux
flux_4 = cal.flux(10**log_E_e, 10**log_E_b, p, 10**log_dNe, 10**log_E_iso_n, 10**log_Eta_0_n, z, theta_j_n, theta_v_fj, t1)
flux_3 = cal.flux(10**log_E_e, 10**log_E_b, p, 10**log_dNe, 10**log_E_iso_w, 10**log_Eta_0_w, z, theta_j_w, theta_v_fj, t1)
flux_2 = cal.flux(10**log_E_e, 10**log_E_b, p, 10**log_dNe, 10**log_E_iso_n, 10**log_Eta_0_n, z, theta_j_n, theta_v_rj, t1)
flux_1j = cal.flux(10**log_E_e, 10**log_E_b, p, 10**log_dNe, 10**log_E_iso_w, 10**log_Eta_0_w, z, theta_j_w, theta_v_rj, t1)
# 汇总flux
flux_olap = np.transpose(flux_4 + flux_3 + flux_2 + flux_1)
flux = np.zeros(len(t),dtype=np.float64)
for i in range(len(t)):
mask1 = i < num_r
mask2 = i < num_r + num_o
for j in range(len(t1)):
if t[i] == t1[j]:
flux[i] = np.where(mask1, flux_olap[j, 1], np.where(mask2, flux_olap[j, 2], flux_olap[j, 3]))
return flux
# lnlike函数
def lnlike(theta):
log_E_e, log_E_b, p, log_dNe, log_E_iso_w, log_E_iso_n, log_Eta_0_w, log_Eta_0_n, theta_j_w, theta_j_n, theta_v_fj, theta_v_rj = theta
model = flux_Model(t_all, log_E_e, log_E_b, p, log_dNe, log_E_iso_w, log_E_iso_n, log_Eta_0_w, log_Eta_0_n, theta_j_w, theta_j_n, theta_v_fj, theta_v_rj)
sigma_h = yerr_h
sigma_l = yerr_l
mask = model > y
sigma = np.where(mask, sigma_h, sigma_l)
sqr1 = (y - model) ** 2
sqr2 = sigma ** 2 + model ** 2
num1 = (2 * np.pi * (sigma_l ** 2 + model ** 2)) ** 0.5
num2 = (2 * np.pi * (sigma_h ** 2 + model ** 2)) ** 0.5
return np.sum(- 0.5 * sqr1 / sqr2 - np.log(num1 + num2))
# lnprior函数
def lnprior(theta):
log_E_e, log_E_b, p, log_dNe, log_E_iso_w, log_E_iso_n, log_Eta_0_w, log_Eta_0_n, theta_j_w, theta_j_n, theta_v_fj, theta_v_rj = theta
if (-5.0 < log_E_e < 0.0 and -5.0 < log_E_b < 0.0 and 2.0 < p < 5.0 and -10.0 < log_dNe < 0.0 and 47.0 < log_E_iso_w < 55.0
and 47.0 < log_E_iso_n < 55.0 and 1.0 < log_Eta_0_w < 5.0 and 1.0 < log_Eta_0_n < 5.0 and 0.01 < theta_j_w < np.pi / 2.0
and 0.01 < theta_j_n < np.pi / 2.0 and 0.0 < theta_v_fj < 1.0 and (theta_v_rj + theta_v_fj) == np.pi and theta_j_n < theta_j_w):
return 0.0
return -np.inf
# lnprob函数
def lnprob(theta):
lp = lnprior(theta)
if not np.isfinite(lp):
return -np.inf
if not np.isfinite(lnprior(theta)):
return -np.inf
return lp + lnlike(theta)
# 查找emcee文档
def main():
# 设置nwalkers,以及参数的初始值
pool = Pool(processes=16) # 可修改
ndim, nwalkers, nsteps = 12, 50, 1000
print('Running MCMC...')
initial = [(np.random.uniform(-0.869, -0.868), # log_E_e
np.random.uniform(-2.60, -2.55), # log_E_b
np.random.uniform(2.134, 2.135), # p
np.random.uniform(-1.59, -1.57), # log_dNe
np.random.uniform(50.82, 50.83), # log_E_iso_w
np.random.uniform(51.94, 51.95), # log_E_iso_n
np.random.uniform(2.43, 2.45), # log_Eta_0_w
np.random.uniform(2.95, 2.96), # log_Eta_0_n
np.random.uniform(0.0817, 0.0823), # theta_j_w
np.random.uniform(0.029, 0.031), # theta_j_n
np.random.uniform(0.519, 0.521), # theta_v_fj
np.random.uniform(2.60, 2.62), # theta_v_rj
) for i in range(nwalkers)]
# set up the backend
task = 'backend_task'
filename = './%s/backend.h5' % task
backend = emcee.backends.HDFBackend(filename)
rm_backend = False
if os.path.exists(filename) and rm_backend is False:
print('The backend file already exists, start with nestp = (%s/%s)' % (str(backend.iteration), str(nsteps)))
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, pool=pool, backend=backend)
sampler.run_mcmc(initial, nsteps=nsteps - backend.iteration, progress=True)
else:
if os.path.exists(filename):
os.remove(filename)
print('Start with a new backend file.')
backend.reset(nwalkers, ndim)
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, pool=pool, backend=backend)
sampler.run_mcmc(initial, nsteps, progress=True)
burnin = int(0.4 * (nwalkers * nsteps))
sampler = sampler.flatchain[burnin:, ]
np.savetxt('sampler.txt', sampler, fmt='%s')
log_E_e, log_E_b, p, log_dNe, log_E_iso_w, log_E_iso_n, log_Eta_0_w, log_Eta_0_n, theta_j_w, theta_j_n, theta_v_fj, theta_v_rj = map(lambda v: (v[1], v[2] - v[1], v[1] - v[0]),
zip(*np.percentile(sampler, [15, 50, 85], axis=0)))
result = [log_E_e, log_E_b, p, log_dNe, log_E_iso_w, log_E_iso_n, log_Eta_0_w, log_Eta_0_n, theta_j_w, theta_j_n, theta_v_fj, theta_v_rj]
np.savetxt('parameter_result.txt', result, fmt='%s')
print("The results have been saved!")
```
