Передача числового массива против скаляра в функцию дает противоречивые результаты

Question

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

import numpy as np

def fn(n, x0):
    return (np.sin(np.arcsin(x0**0.5)*2**n))**2

n = np.arange(100)

x0 = 0.3
print(fn(n, x0))
print(fn(50, x0))

Сгенерированный результат:

[0.3        0.84       0.5376     0.99434496 0.02249224 0.08794536
 0.32084391 0.87161238 0.44761695 0.98902407 0.04342185 0.16614558
 0.55416492 0.98826465 0.04639054 0.17695382 0.58256466 0.97273231
 0.10609667 0.37936067 0.94178461 0.21930545 0.68484228 0.86333333
 0.47195556 0.99685404 0.01254426 0.04954761 0.18837059 0.61154843
 0.95022779 0.18917976 0.         0.         0.         0.
 0.         0.         0.         0.         0.         0.
 0.         0.         0.         0.         0.         0.
 0.         0.         0.         0.         0.         0.
 0.         0.         0.         0.         0.         0.
 0.         0.         0.         0.         0.         0.
 0.         0.         0.         0.         0.         0.
 0.         0.         0.         0.         0.         0.
 0.         0.         0.         0.         0.         0.
 0.         0.         0.         0.         0.         0.
 0.         0.         0.         0.         0.         0.
 0.         0.         0.         0.        ]
0.9931071163166798

Два результата несовместимы, потому что элемент индекса 50 (где n =50) в результатах fn(n, x0) равен нулю, тогда как оценка функции для 50 приводит к ненулевому значению.Почему это так?

Почему большинство элементов результирующего массива имеют f(n, x0) нули?Это не должно быть в соответствии с математикой.

Answer 1

Это требует некоторой отладки, которая довольно типична для работы с numpy.Пара общих причин:

1. Целочисленные типы данных, страдающие от переполнения или целочисленного деления
2. Типы данных с плавающей запятой, дающие неточное поведение
3. Вычисления, вызывающие nan или inf перейти к последующим шагам

Способ выяснить, какая причина (или если это совершенно другая причина) состоит в том, чтобы пройти шаг за шагом.

>>> x0 = 0.3
>>> x0**0.5
0.5477225575051661

Здесь нет проблем,Не удивительно, так как этот шаг не связан с n.Следующий шаг

 >>> 2**n
array([                   1,                    2,                    4,
                          8,                   16,                   32,
                         64,                  128,                  256,
                        512,                 1024,                 2048,
                       4096,                 8192,                16384,
                      32768,                65536,               131072,
                     262144,               524288,              1048576,
                    2097152,              4194304,              8388608,
                   16777216,             33554432,             67108864,
                  134217728,            268435456,            536870912,
                 1073741824,           2147483648,           4294967296,
                 8589934592,          17179869184,          34359738368,
                68719476736,         137438953472,         274877906944,
               549755813888,        1099511627776,        2199023255552,
              4398046511104,        8796093022208,       17592186044416,
             35184372088832,       70368744177664,      140737488355328,
            281474976710656,      562949953421312,     1125899906842624,
           2251799813685248,     4503599627370496,     9007199254740992,
          18014398509481984,    36028797018963968,    72057594037927936,
         144115188075855872,   288230376151711744,   576460752303423488,
        1152921504606846976,  2305843009213693952,  4611686018427387904,
       -9223372036854775808,                    0,                    0,
                          0,                    0,                    0,
                          0,                    0,                    0,
                          0,                    0,                    0,
                          0,                    0,                    0,
                          0,                    0,                    0,
                          0,                    0,                    0,
                          0,                    0,                    0,
                          0,                    0,                    0,
                          0,                    0,                    0,
                          0,                    0,                    0,
                          0,                    0,                    0,
                          0])

Четко показывает проблемное поведение.Там есть куча нулей, перед которыми стоит большое отрицательное число.Если присмотреться более внимательно к результату

>>> (2**n).dtype
dtype('int64')

, вы можете увидеть, что от 2 до n-го числа все еще является целым числом.Похоже, что происходит простой сдвиг битов, приводящий ко всем нулям, когда результат больше не может содержаться в целочисленном типе данных.

Исправление простое, используйте n.astype(float) вместо n.

Answer 2

Пожалуйста, имейте в виду, что простая замена целых чисел числами с плавающей точкой здесь не является правильным решением.На самом деле значения, возвращаемые этим методом для больших n, являются полностью поддельными.

Позвольте мне продемонстрировать это, используя bigfloat библиотеку, обеспечивающую произвольные значения с плавающей точкой.

Используя float64, возвращаемое значениедля n = 99 x0 = 0,3

>>> fn(99., 0.3)
0.013782590413701074

равно ~ 0,0138.

Теперь мы вычислим правильное значение, используя невероятно высокую точность 1000 бит.

>>> def fnbig(n, x0):
...     n, x0 = map(bigfloat.BigFloat, (n, x0))
...     return bigfloat.pow(bigfloat.sin(bigfloat.asin(bigfloat.sqrt(x0))*bigfloat.pow(2, n)), 2)
... 
>>> 
>>> bigfloat.setcontext(bigfloat.Context(1000))

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

>>> fnbig('99', '0.3')
BigFloat.exact('0.363780859940401348053691101648398065131477584225708461696799538248050278540782181716110363889498612214432889606382752875154011855764448898240841915231368492158238806206980341185053867226372528105024157964509865633147960964164133657255856469376571664623973084231004713906743471127849494395877727320492003', precision=1000)

Возвращает ~ 0.364.

Чтобы убедиться, что это правильно, давайте удвоим точность до 2000 бит.

>>> bigfloat.setcontext(bigfloat.Context(2000))
>>> 
>>> fnbig('99', '0.3')
BigFloat.exact('0.3637808599404013480536911016483980651314775842257084616967995382480502785407821817161103638894986122144328896063827528751540118557644488982408419152313684921582388062069803411850538672263725281050241579645098656331479609641641336572558564693765716646239730842310047139067975945820088206827783571320258628882284629795545097600685961974610320482001970915733612836861863674071009032317962504679512051859460424746278292327826581975723619660002116915303723311156451829258099225827808017028059470793304713100650332080089174169114171398280313842625628566029927379227478504732491009738418661061753082431884081337', precision=2000)

Значения, возвращаемые с точностью 1000 и 2000, по сути, одинаковы, поэтому мы можем быть уверены, что ониправильно.

Значение, возвращаемое с использованием "нормальной" арифметики с плавающей точкой, напротив, в значительной степени случайное число.

Answer 3

Для вычисления математических выражений, в которых в результатах ожидаются десятичные значения, важно специально определить числовые значения в функции как числа с плавающей запятой:

def fn(n, x0):      
  return (np.sin(np.arcsin(x0**0.5)*2.**n))**2.

Затем для (fn(n, x0)) я получаюрезультат:

[0.3        0.84       0.5376     0.99434496 0.02249224 0.08794536
 0.32084391 0.87161238 0.44761695 0.98902407 0.04342185 0.16614558
 0.55416492 0.98826465 0.04639054 0.17695382 0.58256466 0.97273231
 0.10609667 0.37936067 0.94178461 0.21930545 0.68484228 0.86333333
 0.47195556 0.99685404 0.01254426 0.04954761 0.18837059 0.61154843
 0.95022779 0.18917976 0.61356311 0.94841368 0.19570069 0.62960771
 0.93280737 0.25071114 0.75142025 0.74715144 0.75566467 0.7385423
 0.77239028 0.70321414 0.83481605 0.55159286 0.98935271 0.04213571
 0.16144118 0.5415117  0.99310712 0.02738149 0.10652697 0.38071589
 0.9430852  0.21470202 0.67442025 0.87831031 0.42752525 0.97898964
 0.08227569 0.3020256  0.84322454 0.52878765 0.99668508 0.01321571
 0.05216421 0.19777242 0.63463397 0.92749478 0.26899286 0.78654281
 0.67157288 0.88225099 0.41553673 0.97146383 0.11088744 0.39436567
 0.95536555 0.17056885 0.56590047 0.98262851 0.06827888 0.25446749
 0.75885515 0.73197605 0.78474845 0.67567327 0.8765556  0.43282351
 0.98194928 0.07089957 0.26349128 0.7762545  0.69473381 0.84831498
 0.51470671 0.99913485 0.0034576  0.01378259]

И (fn(n, x0))[50] оценивается как 0.993107116317.

Хотя я получаю эти результаты с n = np.arange(100), хорошей практикой будет сохранять все как плавающие и использовать n = np.arange(100.).