Как создать график доверительного контура для интегральной функции фитинга Хи-квадрат с различными параметрами?

Question

Я пытаюсь воспроизвести график, похожий на этот (http://supernova.lbl.gov/PDFs/SCP2003OmegaConfTrans.pdf), показывая доверительные области для значений Omega_M и Omega_Lambda.

У меня есть три функции: первая возвращаетинтегрировать для второго, который выполняет интеграцию по scipy. Оба принимают аргумент z, и подинтегральная функция содержит оба Omegas в выражении, но они не принимаются в качестве аргументов. Третья функция возвращает соответствие квадрата Хи для созданных величин моделииз данных, сравнивая их с измеренными величинами.

Я пытаюсь создать двумерный массив значения хи-квадрат модели для каждой комбинации Omega Matter и Omega Lambda в заданном диапазоне (0,0-3,0 и -1,0-3,0 соответственно), а затем нанесите его на график как контурное изображение.

Я попытался использовать двойной цикл for, повторяя длины двух используемых пространств np.linspaceдля генерации диапазона значений Omega, затем установите для Omegas значения ith и jth в lists, вызывая функцию квадрата Хи и устанавливая индекс [i, j] в пустом массиве, равный значению, возвращаемому функцией.

Затем я нанес это с помощью plt.contourf против np.meshgrid, созданного из двух np.linspaces пробных значений.Кажется, что это работает, но возвращает несколько значений nan и не дает сюжет, очень похожий на тот, который я связал.

Omega_Matter_Trials = np.linspace(0.0, 3.0, 100)
Omega_Lambda_Trials = np.linspace(-1.0, 3.0, 100)
Omega_Matter = Omega_Matter_Best
Omega_Lambda = Omega_Lambda_Best

x, y = np.meshgrid(Omega_Matter_Trials, Omega_Lambda_Trials)


# Function to be integrated to find best R0eta using best Om and Ol
def H2(z):
    return (H0_SI*(Omega_Matter*z*(1+z)**2.0 + Omega_Lambda*(1- (1+z)**2.0)         
    +(1+z)**2.0)**0.5)**-1.0
# Function returning integral equal to R0eta
def hr_R0eta(z):
    integral = integrate.quad(H2, 0.0, z)
    #print integral
    return c*integral[0]

def getChi2():
    Chi2 = [] 
    for r in range(len(all_redshifts)):
       model_f = Lpeakbest / (4.0 * np.pi * hr_R0eta(all_redshifts[r])**2.0 
       * (1.0 + all_redshifts[r])**2.0)
       model_m = m0 - 2.5*np.log10(model_f)
       Chi2val = (((all_mags[r]-model_m)/ 
       all_errors[r])**2.0)/len(all_mags)
       Chi2.append(Chi2val)
return np.sum(Chi2)

height = np.zeros_like(x)
for i in range(len(Omega_Matter_Trials)):
    for j in range(len(Omega_Lambda_Trials)):
        Omega_Matter = Omega_Matter_Trials[i]
        Omega_Lambda = Omega_Lambda_Trials[j]
        #print Omega_Matter
        #print Omega_Lambda
        height[i,j] = getChi2() 

levels = np.linspace(0,10)

plt.figure()
plt.contourf(x, y, height, levels)
plt.xlabel('Omega Matter')
plt.ylabel('Omega Lambda')
plt.show()

Я ожидаю, что что-то похожее на этот график http://supernova.lbl.gov/PDFs/SCP2003OmegaConfTrans.pdf,, так как мой код использует те же данные, а функции были использованы в других местах для генерации хорошо подходящей модели с использованием значений наилучшего соответствия обоих Omegas.Я также ожидаю, что минимум графика будет в точке этих наиболее подходящих значений Omega_Lambda = 0,713 и Omega_M = 0,287.В настоящее время я получаю сюжет , который не имеет минимума в нужном месте и имеет форму, отличную от указанной выше.Если есть способ улучшить это так, чтобы он отображал уровни, относящиеся к доверительным интервалам, и имел правильный минимум, то я был бы очень признателен за это.