# -*- coding: utf-8 -*-

"""
@Datetime: 2019/3/29
@Author: Zhang Yafei
"""
# http://itfish.net/article/61104.html

from functools import reduce

import numpy as np
import sympy
from scipy.stats import bernoulli, binom


def bernoulli_mle():
    # 生成样本
    p_1 = 1.0 / 2  # 假设样本服从p为1/2的bernouli分布
    fp = bernoulli(p_1)  # 产生伯努利随机变量
    xs = fp.rvs(100)  # 产生100个样本
    print(xs[:30])  # 看看前面30个
    # [0 1 1 1 1 0 0 0 0 1 0 1 0 0 0 0 1 1 1 1 0 1 1 0 0 1 0 0 0 1]

    # 估计似然函数
    x, p, z = sympy.symbols('x p z', positive=True)
    phi = p ** x * (1 - p) ** (1 - x)  # 分布函数
    L = np.prod([phi.subs(x, i) for i in xs])  # 似然函数
    print(L)
    # p**51*(-p + 1)**49

    logL = sympy.expand_log(sympy.log(L))
    sol = sympy.solve(sympy.diff(logL, p), p)
    print(sol)
    # 51/100

    b = binom(100, .5)
    # 以均值为中心对称的加总概率
    g = lambda x: b.pmf(np.arange(-x, x) + 50).sum()
    print(g(10))

    b.pmf(np.arange(40, 60))


"""
https://blog.csdn.net/buptgshengod/article/details/38817621
考虑一个抛硬币的例子。假设这个硬币正面跟反面轻重不同。我们把这个硬币抛80次（即，
我们获取一个采样并把正面的次数记下来，正面记为H，反面记为T）。并把抛出一个正面的概率记为，
抛出一个反面的概率记为（因此，这里的即相当于上边的）。假设我们抛出了49个正面，31个反面，
即49次H，31次T。假设这个硬币是我们从一个装了三个硬币的盒子里头取出的。这三个硬币抛出正面的概率
分别为, , .这些硬币没有标记，所以我们无法知道哪个是哪个。使用最大似然估计，通过这些试验数据
（即采样数据），我们可以计算出哪个硬币的可能性最大。这个似然函数取以下三个值中的一个：

我们可以看到当p=2/3时，似然函数取得最大值。这就是p的最大似然估计。

其实就是把实验结果的出现概率写作p的一个函数，然后求该函数较大的p值
"""


def Factorial(x):
    """
    阶乘
    :param x:
    :return:
    """
    return reduce(lambda x, y: x * y, range(1, x + 1))


def maximum_like_lihood(p, sample, x):
    """ 极大似然估计 """
    f1 = Factorial(sample) / (Factorial(x) * Factorial(sample - x))
    f2 = (p ** x) * ((1.0 - p) ** (sample - x))

    return f1 * f2


if __name__ == '__main__':
    mle = maximum_like_lihood(p=2.0/3, sample=60, x=40)
    print(mle)