Статистика в Python - PullRequest
Новая
       60

Статистика в Python

0 голосов
/ 25 октября 2018

Я беру урок машинного обучения, и нам дают наше первое статистическое упражнение - программирование.

Итак, упражнение выглядит следующим образом:

Вспомните историю из лекции «Два продавца в Amazon имеют одинаковую цену.Один имеет 90 положительных и 10 отрицательных отзывов.Другой 2 положительный и 0 отрицательный.У кого покупать? »Запишите последующие вероятности надежности (как в лекции).Вычислить p (x> y | D1, D2), используя численное интегрирование.Вы можете генерировать бета-версии распределенных сэмплов с помощью функции scipy.stats.beta.rvs (a, b, size) .

. Из лекции известно следующее:

применил две бета-биномиальные модели: p (x | D1) = бета (x | 91, 11) и p (y | D2) = бета (y | 3, 1)

вычислитьвероятность того, что продавец 1 более надежен, чем продавец 2: 

p(x > y | D1, D2 ) = ∫∫ [x > y] Beta (x| 91, 11) Beta (y| 3, 1) dx dy

Так что мои попытки в Python таковы: 

In [1]: import numpy as np
        from scipy import integrate, stats
In [2]: f = lambda x, y: stats.beta.rvs(91, 11, x) * stats.beta.rvs(3, 1, y)
In [3]: stats.probplot(result, x > y)

И я получаю сообщение об ошибке:

... The maximum number of subdivisions (50) has been achieved....

но в конечном итоге есть ответ на расчет, который составляет ок.1.7.(Нам говорят, что ответ составляет приблизительно 0,7)

Мой вопрос: Как рассчитать часть [x> y], что означает: вероятность того, что продавец 1 (x)надежнее продавца 2 (у)?

1 Ответ

0 голосов
/ 26 октября 2018

почти верно, я бы сделал что-то вроде: 

from scipy import stats

N = 10000
P = sum(stats.beta.rvs(3, 1, size=N) < stats.beta.rvs(91, 11, size=N))
P / N

, и если вы хотите графическое отображение: 

import matplotlib.pyplot as plt
import numpy as np

X = np.linspace(0.6, 0.8, 501)
Y = stats.beta.pdf(X, 1 + P, 1 + N - P)

plt.plot(X, Y)

может быть библиотечный код для более удобного построения графика.

Выше приведена оценка ответа по методу Монте-Карло.Если вам нужна более точная числовая оценка, вы можете получить ее в квадратуре, используя: 

from scipy.integrate import dblquad
from scipy import stats

a = stats.beta(91, 11)
b = stats.beta(3, 1)

dblquad(
    lambda x, y: a.pdf(x) * b.pdf(y),
    0, 1, lambda x: x, lambda x: 1)

, что дает мне оценку ~ 0,712592804 (с ошибкой 2e-8).

Если выхочу стать еще точнее, вы захотите сделать что-то аналитическое:

https://stats.stackexchange.com/questions/7061/binomial-probability-question

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

...