Уточнить математическую формулу в коде

У меня есть функция для вычисления нормального распределения в Python:

def norm_cdf(z):
  """ Use the norm distribution functions as of Gale-Church (1993) srcfile. """
  # Equation 26.2.17 from Abramowitz and Stegun (1964:p.932)

  t = 1.0 / (1+0.2316419*z) # t = 1/(1+pz) , p=0.2316419
  probdist = 1 - 0.3989423*math.exp(-z*z/2) * ((0.319381530 * t)+ \
                                         (-0.356563782* math.pow(t,2))+ \
                                         (1.781477937 * math.pow(t,3)) + \
                                         (-1.821255978* math.pow(t,4)) + \
                                         (1.330274429 * math.pow(t,5)))
  return probdist

Но я должен придерживаться поля PEP8 и 80, следовательно, уродливый \ s. Как еще я должен префикс моего кода?

В математической форме

$$ \ {начать выравнивать *} \ textrm {norm_cdf} (z) = 1 - 0,3989423 e ^ \ frac {-z ^ 2} {2} (& 1,330274429 t ^ 5 - 1,821255978 t ^ 4 \\ & amp ;: 1,781477937 t ^ 3 - 0,356563782 t ^ 2 + 0,319381530 t) \ Конец {выравнивать *} $$

, где

$$ t = \ frac {1} {1 + 0.2316419 z} $$

42 голоса | спросил alvas 7 PM00000050000002931 2014, 17:12:29

8 ответов


40

Позвольте мне процитировать замечательную книгу Numerical Recipes в C ++ (но также применимую):

  

Мы предполагаем, что вы знаете достаточно, чтобы никогда не оценивать полином таким образом:

p=c[0]+c[1]*x+c[2]*x*x+c[3]*x*x*x+c[4]*x*x*x*x;
     

или (еще хуже!),

p=c[0]+c[1]*x+c[2]*pow(x,2.0)+c[3]*pow(x,3.0)+c[4]*pow(x,4.0);
     

Придет (компьютерная) революция, все лица, признанные виновными в таком преступном поведении, будут выполняться в кратчайшие сроки, а их программ не будет!

(Вы можете найти страницу в своем издании в аналитическом индексе под кадром ", особенно плохо" . Мне нравится эта книга.)

Есть две причины не делать этого: точность и производительность. Правильный способ вычисления многочлена таков:

-t * (0.319381530  +  t * (-0.356563782 + t * (1.781477937 + t * (-1.821255978 + 1.330274429 * t))))

И вы можете, конечно, расколоться по вашему усмотрению, так как строки новой строки в скобках игнорируются. Помните PEP: «Предпочтительное место для разрыва вокруг двоичного оператора - после оператора, а не до него».

-t * (0.319381530  +  t * (-0.356563782 +
    t * (1.781477937 + t * (-1.821255978 + 1.330274429 * t))))

Другой альтернативой является сохранение коэффициентов в списке:

coeff = [0, 0.319381530, -0.356563782, 1.781477937, -1.821255978, 1.330274429]
poly = coeff[-1]
for c in coeff[-2::-1]:
    poly *= x
    poly += c

Я выполняю операции, чтобы избежать выделения и освобождения памяти, но это имеет смысл только в том случае, если x является массивом NumPy. Если вы оцениваете один номер, вы можете просто использовать более приятное выражение:

poly = poly * x + coeff[i]

Но я бы придерживался первого, потому что он более общий.

Конечно, результат должен быть умножен на префактор:

return 1 - 0.3989423*math.exp(-z*z/2) * poly

Или, если вы хотите сделать это на месте:

z2 = z * z # Be careful not to modify your input!
z2 *= 0.5  # Multiplication is faster than division.
np.exp(z2, out=z2)

probd = z2 * poly
probd *= -0.3989423
probd += 1
return probd

Бонус-трек:

Если вы применяете эту функцию к большим массивам (более тысячи номеров), вы можете воспользоваться первой техникой в ​​numexpr:

expr += '1 - 0.3989423* exp(-z * z / 2) * '
expr += '(-t * (0.319381530  +  t * (-0.356563782 +  t * '
expr += '(1.781477937 + t * (-1.821255978 + 1.330274429 * t)))))'
ne.evaluate(expr)

Это будет компилировать выражение для вас и прозрачно использовать столько ядер, сколько у вас есть.

ответил Davidmh 7 PM00000070000003731 2014, 19:55:37
22

Как выясняется, вопрос о аналогичном вопросе был задан недавно на Math.SE. Вместо , используйте встроенные функции в Python.

Ваш norm_cdf(z) является просто численным приближением для

$$ P (z) = \ frac {1} {\ sqrt {2 \ pi}} \ int _ {- \ infty} ^ {z} e ^ {- t ^ 2/2} \ dt = \ int _ {- \ infty} ^ {z} Z (t) \ dt = \ frac {1} {2} \ left (1 + \ mathrm {erf} \ left (\ frac {z} {\ sqrt {2}} \ right) \ right) = \ frac {1} {2} \ mathrm {erfc} \ left (- \, \ frac {z} {\ sqrt {2}} \ right) $$

Поэтому вы можете просто использовать math.erfc() (доступно с Python 2.7) и получить более точный результат (особенно для очень отрицательных значений \ $ z \ $).

import math

def norm_cdf(z):
    return 0.5 * math.erfc(-x / math.sqrt(2))

Еще лучше, просто используйте scipy.stats.norm.cdf() !

ответил 200_success 7 PM00000080000002231 2014, 20:56:22
18

Я буду рассматривать только так называемые «магические числа», о которых упомянули несколько рецензентов.

Иногда, когда вы работаете в чистой математике, то, что на первый взгляд кажется «магическим числом», действительно не . Может быть, сами цифры являются лишь частью заявления о проблеме. Я думаю, этот вопрос сводится к следующему: вы можете придумать имя, которое более наглядно, чем число? Если есть хорошее имя, вы, вероятно, должны его использовать.

На первый взгляд, я думал, что ваши номера являются неотъемлемой частью проблемы. Но когда я посмотрел на Абрамовица и Стегуна, я увидел, что ссылочная формула уже назвала ваши уродливые константы. Имена: p (которые вы упомянули в комментарии), и b1 через b5. Вы должны использовать эти имена в коде, потому что они создают очень четкую ссылку на исходное определение формулы.

Когда вы решили, что было бы неплохо добавить комментарий p=0.2316419, это было очень убедительное доказательство того, что имя должно . (И как только code говорит p=0.2316419, комментарий должен быть удален.)

Кстати, было очень хорошо, чтобы вы включили точную ссылку Abramowitz и Stegun в комментарий.

ответил GraniteRobert 8 AM000000120000000831 2014, 00:40:08
16

Вместо math.pow используйте встроенный ** оператор. Вы не нуждаетесь в \ s в EOL, потому что круглые скобки, окружающие выражение, позволяют ему неявно охватывать несколько строк. Поэтому после внесения обоих этих изменений я завершаю:

def norm_cdf(z):
    """ Use the norm distribution functions as of Gale-Church (1993) srcfile. """
    # Equation 26.2.17 from Abramowitz and Stegun (1964:p.932)

    t = 1.0 / (1 + 0.2316419*z) # t = 1/(1+pz) , p=0.2316419
    probdist = 1.0 - (   (0.319381530  * t)
                       + (-0.356563782 * t**2)
                       + (1.781477937  * t**3)
                       + (-1.821255978 * t**4)
                       + (1.330274429  * t**5)) * 0.3989423*math.exp(-z*z/2)
    return probdist

Я также перегруппировал одно из умножений, чтобы сделать приоритет более очевидным и читаемым.

После этих изменений единственная проблема, которую я все еще вижу, - это все магические числа. Я не знаю, как эти константы прибывают, но это может помочь читаемости, если константам могут быть предоставлены значимые имена. Иногда с формулами не очень много значимого имени.

ответил godlygeek 7 PM00000060000005531 2014, 18:56:55
15

Не уверен, что это помогает, но вы можете легко определить функцию для оценки значения полинома в данной позиции

def evaluate_polynom(pol, x):
    return sum(a * math.pow(x, i) for i, a in enumerate(pol))

Тогда

(0.319381530 * t) + (-0.356563782* math.pow(t,2)) + (1.781477937 * math.pow(t,3)) + (-1.821255978* math.pow(t,4)) + (1.330274429 * math.pow(t,5))

становится:

evaluate_polynom([0, 0.319381530, -0.356563782, 1.781477937, -1.821255978, 1.330274429], t)
ответил Josay 7 PM00000050000004331 2014, 17:33:43
6

В случае опасности (дальнейшего) раздражения модераторов я посоветую вам следовать советам о численных методах, которые исходят из любой из книг Numeric Recipes 1 .

Хотя метод Хорнера может помочь в некоторых случаях, это долгий путь от панацеи. Вместо того, чтобы фиксировать проблему округления при суммировании значений, он пытается избежать этой проблемы. К сожалению, это лишь частично успешно. Для большинства полиномов все равно будут входные данные, для которых результаты будут относительно бедными, даже в лучшем случае.

Если значения, выражаемые членами полинома, могут приводить к числовой нестабильности при суммировании, вы можете рассмотреть возможность генерации каждого отдельно, затем используя что-то вроде Kahan sum , чтобы суммировать эти условия. Если вам это интересно, это также может дать вам поле погрешности вместе с самой суммой.

Вероятно, лучше еще использовать Langlois, et al. компенсировал схему Хорнера . По крайней мере, в последний раз, когда я смотрел внимательно, это, по-видимому, было в значительной степени актуальным для оценки полиномов. Он поддерживает примерно ту же точность результата, что и вы, используя схему Хорнера с использованием чисел с плавающей запятой, с двойной точностью (например, с использованием 64-битного двойника, он дает примерно такую ​​же точность, как схема Хорнера с 128-битным квад- точность плавающей запятой, но без ограничения скорости, которое обычно несут). Подобно сумме Kahan, это поддерживает вычисление ошибки, связанной с результатом.


1. Для справки см. Критические замечания, такие как: http : //www.stat.uchicago.edu/~lekheng/courses/302/wnnr/nr.html
http://www.lysator.liu.se/c/num-recipes-in-c.html
Я думаю, что более точное резюме, чем «замечательно», следующее: «Я обнаружил, что Numerical Recipes предоставляет достаточно информации для человека, чтобы попасть в неприятности, потому что после прочтения NR каждый думает, что кто-то понимает что происходит." Суб>

ответил Jerry Coffin 10 PM00000010000001631 2014, 13:40:16
2

Я буду честным здесь: я полностью сосать на питоне.

Тем не менее, можно ли объявить некоторые из этих числовых литералов константами? Это очистит ваш код и еще раз прояснит сам код.

ответил Pimgd 7 PM00000050000001931 2014, 17:21:19
2

тогда:

a*x**3 + b*x**2 + c*x = ((a*x + b)*x + c)*x

h /t @ Эмили Л. для справки «Хорнер»:

https://en.wikipedia.org/wiki/Horner%27s_method

и h /t @ Davidmh за то, что заметили улучшения в скорости вычислений /точности этого метода.

gale-church процитировала это так в 1990 году:

import math

def pnorm(z):

    t = 1 / (1 + 0.2316419 * z)
    pd = (1 - 0.3989423 *
      math.exp(-z * z / 2) *
        ((((1.330274429 * t - 1.821255978) * t
           + 1.781477937) * t - 0.356563782) * t + 0.319381530) * t)

    return pd

Этот метод удобно избегает проблемы t ^ n.

цитата:

 введите описание изображения здесь>> </a> </p>

<p> <a href= введите описание изображения здесь>> </a> </p>

<p> Источник: </p>

<p> <a href= http://www.aclweb.org/anthology/J93-1004

стр. 21 из 28 в формате pdf

стр. 95 журнала Вычислительная лингвистика Том 19, номер 1

Я могу «прикрыть» до:

def pnorm(z):

t = 1 / (1 + 0.2316419 * z)
pd = (1 - 0.3989423 * math.exp(-z * z / 2) *
      ((((1.330274429 * t - 
          1.821255978) * t + 
          1.781477937) * t - 
          0.356563782) * t + 
          0.319381530) * t )

return pd

, если вы проверите

Abromowitz and Stegun, Руководство по математическим функциям

стр. 932 уравнение 26.2.17

цитата:

http://people.math.sfu.ca/~cbm /aands/page_932.htm

вы найдете следующее:

 введите описание изображения здесь>> </a> </p>

<p>, из которого мы можем создать таблицу, дающую нам: </p>

<pre><code>def pnorm (z):

    p = 0,2316419
    b1 = 0,319381530
    b2 = -0.356563782
    b3 = 1,781477937
    b4 = -1,821255978
    b5 = 1,330274429
    t = 1 /(1 + p * z)
    pd = (1 - 0,3989423 * math.exp (-z * z /2) *
          ((((b5 * t + b4) * t + b3) * t + b2) * t + b1) * t)

    return pd
</code></pre>

<p> Затем с предыдущей страницы; 931 вы найдете: </p>

<p> <a href= введите описание изображения здесь>> </a> </p>

<pre><code>Zx = (1 /âš (2 * Ï €)) * e (-z * z /2)
</code></pre>

<p> в python: </p>

<pre><code>Zx = (1 /math.sqrt (2 * math.pi)) * math.exp (-z * z /2)
</code></pre>

<p> и находим, что (1 /âš (2 * Ï €)) = 0,3989423 </p>

<p> тоже, мне кажется, мне это нравится: </p>

<pre><code>t * (b1 + t * (b2 + t * (b3 + t * (b4 + t * b5))))
</code></pre>

<p> лучше, чем: </p>

<pre><code>(((b5 * t + b4) * t + b3) * t + b2) * t + b1) * t
</code></pre>

<p>, значит, наконец: </p>

<pre><code>import math

def pnorm (z):

    p = 0,2316419
    b1 = 0,319381530
    b2 = -0.356563782
    b3 = 1,781477937
    b4 = -1,821255978
    b5 = 1,330274429
    t = 1 /(1 + p * z)
    Zx = (1 /math.sqrt (2 * math.pi)) * math.exp (-z * z /2)
    pd = Zx * t * (b1 + t * (b2 - t * (b3 + t * (b4 - t * b5))))

    return (1 - pd)
</code></pre>

<p> проверка моей работы против op </p>

<pre><code>import matplotlib.pyplot как plt
импортировать numpy как np
импортная математика




def norm_cdf (z):
  

 введите описание изображения здесь>> </a> </p>

<p> Я ожидал, что метод Хорнера будет быстрее, поэтому я провел тест времени, заменив: </p>

<pre><code>#Zx = (1.0 /math.sqrt (2.0 * math.pi)) * math.exp (-z * z /2.0)
Zx = 0,3989423 * math.exp (-z * z /2.0)
</code></pre>

<p>, чтобы сделать его честным и повысить разрешение np.arrange до0,0001: </p>

<pre><code>t0 = time.time ()
для z в np.arange (-3,3,0.0001):
    zf = pnorm (z)
t1 = время.time ()
для z в np.arange (-3,3,0.0001):
    zf = norm_cdf (z)
t2 = time.time ()

print ('pnorm time:% s'% (t1-t0))
print ('norm_cdf time:% s'% (t2-t1))
</code></pre>

<p>, и результаты, вращающие мой четырехъядерный процессор AMD 7950 FM2 + w /16G, довольно сильно (хотя и с несколькими другими приложениями) ... не оправдали мои ожидания: </p>

<pre><code>> > >
Время pnorm: 81.4725670815
norm_cdf время: 80.7865998745
</code></pre>

<p> Метод Хорнера не был быстрее </p></body></html>

ответил litepresence 20 MaramMon, 20 Mar 2017 05:34:29 +03002017-03-20T05:34:29+03:0005 2017, 05:34:29

Похожие вопросы

Популярные теги

security × 330linux × 316macos × 2827 × 268performance × 244command-line × 241sql-server × 235joomla-3.x × 222java × 189c++ × 186windows × 180cisco × 168bash × 158c# × 142gmail × 139arduino-uno × 139javascript × 134ssh × 133seo × 132mysql × 132