Я в настоящее время следую за Think Bayes, вводным текстом к Байесовской статистике (http://www.greenteapress.com/thinkbayes/thinkbayes.pdf) и пытался повторить главу 10. В этой главе набор данных для высот используется для генерации байесовской модели для населения. среднее значение и стандартное отклонение для конкретного образца.
Идея состояла в том, чтобы сделать чистые вычислительные методы более устойчивыми с помощью различных способов, включая приблизительные байесовские вычисления, вероятности записи в журнал и ограничение диапазона предшествующих значений mu и sigma для результирующего распределения с использованием известных правил выборки (например, S = sigma /sqrt(n)).
Мой код, кажется, должен работать, но он дает мне ответы, которые явно неверны:
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
def MakeRange(estimate, stderr, num_stderrs,num_points=100):
spread = stderr * num_stderrs
array = np.linspace(estimate-spread,estimate+spread,num_points)
return array
def FindPriorRanges(xs, num_stderrs=3.0):
# compute m and s
n = len(xs)
m = np.mean(xs)
s = np.std(xs)
# compute ranges for m and s
stderr_m = s / np.sqrt(n)
mus = MakeRange(m, stderr_m, num_stderrs)
stderr_s = s / np.sqrt(2 * (n-1))
sigmas = MakeRange(s, stderr_s, num_stderrs)
return mus, sigmas
def Likelihood(x,mus,sigmas):
like=np.nan_to_num(np.fromfunction(lambda x,a,b:stats.norm.pdf(x,mus[a],sigmas[b]), (len(x),len(mus),len(sigmas)),dtype=np.int))#stats.norm.pdf(x,mu,sigma)
return like
def logLikelyhood(x,mu,sigma):
like=Likelihood(x, mu, sigma)
largest=np.max(like)
return np.nan_to_num(np.log(like/largest))
#his code also has something to remove 0 probability elements to prevent -Inf when you apply log
def expLogLikelyhood(logdist):
m=np.max(logdist)
return np.exp(logdist-m)
#note that this is extremely slow, and was made to ensure the algorithm is working exactly as expected
def approxBayesianComputation(data,mus,sigmas):
n=len(data)
m=np.mean(data)
s=np.std(data)
#he adds together the likelyhoods (log) of the mean and the sigma occuring given the particular dataset (rather than the exact data we got here (in terms of numbers))
#would be two seperate apply along axes
#muloglike=np.apply_along_axis(func1d, 0, arr)
#very slow
#using log likelyhood (same result as below)
# arr=np.ones([len(mus),len(sigmas)])
# for muind,mu in enumerate(mus):
# for sigmaind,sigma in enumerate(sigmas):
# stderr_m = sigma / np.sqrt(n)
# loglike=np.log(stats.norm.pdf(m, mu, stderr_m))
# stderr_s=sigma/np.sqrt(2 * (n-1))
# loglike+=np.log(stats.norm.pdf(s, sigma, stderr_s))
# arr[muind][sigmaind]=loglike
#
# return np.exp(arr)
#without log likelyhood
arr=np.ones([len(mus),len(sigmas)])
for muind,mu in enumerate(mus):
for sigmaind,sigma in enumerate(sigmas):
stderr_m = sigma / np.sqrt(n)
loglike=stats.norm.pdf(m, mu, stderr_m)
stderr_s=sigma/np.sqrt(2 * (n-1))
loglike*=stats.norm.pdf(s, sigma, stderr_s)
arr[muind][sigmaind]=loglike
return arr
def Normalize(arr):
return arr/np.sum(arr)
#data=[1,2,3,4,5,6,7,8,2,3,4,5,6,7,3,3,4,4,1,1,1,1,1,1,1,1,1,1,1,1,12,2]*70
#data=[1,4,2,5,6,7,8,-3,-4,-5,-1]
#data=[1,1,1,1,1,1,2,2,2,0,0,0]
data=[1,2,3,4,6]
#data=[-2,0,2,3,4,5]
#data=[7,7,7,7,7,7,7,7,7,7,7,7,7.1]
#data=[0,-1,1,-2,2,-3,3,-4,4,-5,5,-6,6]
#data=[0,0,0,0]
#data=[5,5,5,5]
#all likely ranges for mu and sigma
mus, sigmas=FindPriorRanges(data)
#sigmas=np.linspace(.01,10,100)
#mus=np.linspace(0,30,1000)
prior=np.ones([len(mus),len(sigmas)])
#make 2d distribution
multarray=Likelihood(data,mus,sigmas)
posterior=Normalize(np.multiply(prior,np.prod(multarray,axis=0)))
plt.imshow(posterior, cmap='gist_heat', interpolation='nearest',extent=[min(sigmas),max(sigmas),min(mus),max(mus)])#extent=[min(mus),max(mus),min(sigmas),max(sigmas)]
plt.show()
#using log likelyhoods for larger datasets (note that this works on large datasets where the above fails due to discretization error
multarray=logLikelyhood(data,mus,sigmas)
#the product of logs is the log of the sum
posterior=Normalize(np.multiply(prior,expLogLikelyhood(np.sum(multarray,axis=0))))
plt.imshow(posterior, cmap='gist_heat', interpolation='nearest',extent=[min(sigmas),max(sigmas),min(mus),max(mus)])
plt.show()
plt.imshow(approxBayesianComputation(data, mus, sigmas), cmap='gist_heat', interpolation='nearest',extent=[min(sigmas),max(sigmas),min(mus),max(mus)])
plt.show()
То, что делает эта программа, вычисляет диапазон mus, средств и сигм, стандартных отклонений, на основе конкретной выборки данных (названных данных в коде). В настоящее время он генерирует 100 точек по умолчанию для каждого параметра.
np.fromfunction используется для распространения массива, имеющего размеры (len(data),len(mus),len(sigmas)), используя обычную pdf-функцию scipy (в качестве точечной функции). Затем выполняется умножение вдоль оси данных для обновления распределения, которое затем умножается на предыдущее, которое я имею в качестве равномерного распределения. Для logLikelyhood процесс почти идентичен, за исключением того, что он складывается вдоль оси данных (сумма журналов - это журнал продукта), преобразует журнал обратно в стандартные вероятности путем возведения в степень, а затем умножает на предыдущее.
Однако, когда я запускаю код, я замечаю несколько странных вещей:
- Распределение не центрировано в центре рассчитанных вероятных диапазонов для мю и сигмы, как я считаю, должно быть
- По-видимому, распределение имеет правильный перекос для му, а также расширяющийся «хвост кометы» для сигмы при увеличении му. Обычно вы ожидаете, что он будет идеально эллиптическим.
- Максимум mu и sigma в дистрибутиве - это не mu и sigma данных, что наверняка неправильно.
- Для некоторых закомментированных распределений данных данные, похоже, даже не похожи на упомянутую ранее "комету" или эллипс.
И вероятность, и логоподобие дают одинаковые распределения последовательно. Тем не менее, я получаю некоторое деление на 0 ошибок, а для некоторых дистрибутивов с логарифмической вероятностью я получаю nan значений, которые я сгладил с помощью np.nan_to_num(). Однако даже для дистрибутивов, не имеющих ни одной из этих ошибок (которые, кажется, не оказывают существенного влияния на результат). Кроме того, approxBayesianComputation, предназначенный для прозрачности процесса, который происходил (он медленно умножает все с помощью python для циклов), дает ответы с аналогичными свойствами для более быстрых реализаций numpy.