Проблема с нормализацией функции с помощью модуля интеграции scipy

#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. И я постараюсь склеить ваше предложение и комментарий Уоррена выше. Но большое вам спасибо.