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

Question

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

Итак, упражнение выглядит так:

Вспомните историю из лекции "Два продавца на 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 (y)?

Answer 1

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

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

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