Целочисленная оптимизация/максимизация в numpy

Question

Мне нужно estimate размер популяции, найдя значение n, которое максимизирует scipy.misc.comb(n, a)/n**b, где a и b являются константами. n, a и b - целые числа.

Очевидно, что я мог бы иметь петлю в range(SOME_HUGE_NUMBER), вычислить значение для каждого n и вырваться из цикла, как только достигнет перегиба на кривой. Но я задавался вопросом, был ли элегантный способ сделать это с помощью (например) numpy/scipy или какой-нибудь другой элегантный способ сделать это только в чистом Python (например, как целочисленный эквивалент метода Ньютона?)

Answer 1

As поскольку ваш номер n достаточно мал (меньше 1500), я предполагаю, что самый быстрый способ сделать это - фактически попробовать все возможные значения. Вы можете сделать это быстро, используя numpy:

import numpy as np 
import scipy.misc as misc 

nMax = 1000 
a = 77 
b = 100 
n = np.arange(1, nMax+1, dtype=np.float64) 
val = misc.comb(n, a)/n**b 
print("Maximized for n={:d}".format(int(n[val.argmax()]+0.5))) 
# Maximized for n=181 

Это не особенно элегантно, но довольно быстро для этого диапазона n. Проблема в том, что для n>1484 числитель уже может быть слишком большим для хранения в float. Затем этот метод завершится неудачей, так как вы столкнетесь с переполнениями. Но это не только проблема numpy.ndarray, не работающая с целыми числами python. Даже с ними, вы не смогли бы вычислить:

misc.comb(10000, 1000, exact=True)/10000**1001 

как вы хотите иметь результат с плавающей точкой в ​​вашем делении двух чисел больше, чем максимум в питоне float может держать (max_10_exp = 1024 на моей системе см. sys.float_info().). Вы также не могли использовать ваш range. Если вы действительно хотите сделать что-то подобное, вам придется проявлять большую осторожность численно.

Answer 2

У вас по существу есть приятная плавная функция n, которую вы хотите максимизировать. n требуется, чтобы быть интегралом, но мы можем рассматривать функцию вместо того, чтобы быть функцией от реалов. В этом случае максимизирующее интегральное значение n должно быть близко к (рядом) максимизирующему действительному значению.

Мы можем преобразовать comb в реальную функцию с помощью гамма-функции и использовать численные методы оптимизации, чтобы найти максимум. Другой подход заключается в замене факториалов на приближение Стирлинга. Это дает умеренно сложное, но сговорчивое алгебраическое выражение. Это выражение не трудно отличить и установить на нуль, чтобы найти экстремумы.

Я сделал это и получил

n * (b + (n-a) * log((n-a)/n)) = a * b - a/2 

Это не просто решить алгебраически, но достаточно легко численно (например, с помощью метода Ньютона, как вы предлагаете).

Возможно, я ошибся в алгебре, но я набрал пример a = 77, b = 100 в Wolfram Alpha и получил 180,58, поэтому подход, похоже, работает.