#python #scipy #scipy.stats
Вопрос:
Я хочу нормализовать функцию (например, chi2.pdf scipy) в диапазоне от A до B. Например, chi2.pdf нормализуется для диапазона от 0 до бесконечности, и его интеграл по этой области равен 1. Для этого я могу вычислить интеграл функции от A до B и разделить функцию на этот интеграл. Я могу реализовать это с помощью следующего кода:
import numpy as np from scipy.stats import chi2 from scipy.integrate import quad A = 2 B = 4 df = 3 z = quad(chi2.pdf,A,B,args=(df,A)[0]
Quad передает аргументы df как степень свободы, а A как loc — я хочу, чтобы моя функция хи-квадрат была сдвинута на A по разным причинам. Теперь, когда у меня есть z
, я могу определить новую функцию :
def normalized_chi_2(x,df,A,z): y = chi2.pdf(x,df,A)/z return(y)
Быстрая проверка с помощью интеграции еще раз:
integral_chi2 = quad(normalized_chi_2,A,B,args=(df,A,z)[0] print(integral_chi2) gt;0.9999999999999999
Показывает, что я достиг своей цели. Но наличие двух функций и вычисление Z в основном относительно громоздко, поэтому я решил, что могу определить новую функцию и вычислить Z внутри этой функции.
def normalized_chi_1(x,df,A): z = quad(chi2.pdf,A,B,args=(df,A))[0] y = chi2.pdf(x,df,A) / z return(y)
Теперь, когда я снова выполняю быструю интеграцию:
integral_chi1 = quad(normalized_chi_1,A,B,args=(df,A))[0] print(integral_chi1) gt;0.42759329552910064
Я не получаю 1, и я получаю значение, равное значению исходного ненормализованного chi2.pdf (z выше). Другая проблема заключается в том, что normalized_chi_1
(который принимает df и A и вычисляет свой собственный z) работает очень и очень медленно. Например, метод 2, при котором я вычисляю z вне функции и передаю ее в следующую функцию, занимает ~0,07 секунды, в то время как метод 1, при котором я вычисляю z внутри функции, занимает ~7,30 секунды. В сто раз медленнее.
Комментарии:
1. Вместо интеграции
chi2.pdf
вы могли бы использоватьchi2.cdf(B) - chi2.cdf(A)
. (В конце концов, CDF является неотъемлемой частью PDF.)2. @WarrenWeckesser Это блестяще! Кроме того, поскольку мой chi2.pdf сосредоточен вокруг 2, мне не нужно вычислять chi2.cdf из A, но мне просто нужно chi2.cdf из B. Я протестировал его, и я получаю примерно на 10% быстрее, просто вызывая chi2.cdf(B,df,A), чем вычисление интеграла через quad. Большое спасибо!
Ответ №1:
quad
вероятно, запускает цикл под капотом, и каждый раз, когда ваша функция вызывается, она вызывает другую quad
для вычисления z, все это становится довольно утомительным. Чтобы проверить это, я добавил простую инструкцию печати со счетчиком к исходной функции.
count = 0 def normalized_chi_1(x,df,A): global count z = quad(chi2.pdf,A,B,args=(df,A))[0] print(f"calculating z {count}th time") count = 1 y = chi2.pdf(x,df,A) / z return(y)
Результат, который я получил, был
calculating z 0th time ... calculating z 227th time calculating z 228th time calculating z 229th time calculating z 230th time
Таким образом, вы вычисляете интеграл для z примерно 230
в два раза, что более или менее объясняет 100x
увеличение времени выполнения.
Если вы хотите, чтобы функция вычисляла z, вы можете просто сделать
from functools import lru_cache @lru_cache def get_z(*ars,args): return quad(*ars,args)[0] A = 2 B = 4 df = 3 def normalized_chi_1(x,df,A): z = get_z(chi2.pdf,A,B,args=(df,A)) y = chi2.pdf(x,df,A) / z return(y) integral_chi1 = quad(normalized_chi_1,A,B,args=(df,A))[0]
Это дало мне правильный результат и время работы 0,07 с, но я думаю, что лучше просто определить z в основном.
Комментарии:
1. Я решил, что это зацикливается, но не мог понять, почему он не дает мне правильного ответа. Использование запоминания тоже довольно умно. Я собираюсь сказать, что использование предопределенного дистрибутива scipy имеет преимущество включения cdf. И я постараюсь склеить ваше предложение и комментарий Уоррена выше. Но большое вам спасибо.