dice_example.py

def dice_example_l():
    """
    Example with a biased dice with a previously observed mean
    The lagrange multiplier l is optimized rather than p
    """
    def entropy(l, *args):
        """
        Entropy of p
        """
        p = analytical(l)
        return sum(p*np.log(p))

    def mean_constraint(l, *args):
        """
        Constraint that the weighted mean is fixed
        """
        mean = args[0]
        x = np.arange(1,7)
        p = analytical(l)
        return p.dot(x) - mean

    def analytical(l):
        """
        Convert the lagrange multiplier l to p
        """
        x = np.arange(1,7)
        p = np.exp(-l*x)
        return p/sum(p)

    # Observed mean
    mean = 3
    # Starting guess
    l = -1
    # Set constraints
    constraints = [{'type':'eq',
                    'fun':mean_constraint,
                    'args': (mean,)}
                    ]

    opt = scipy.optimize.minimize(entropy, l, constraints=constraints)
    print(analytical(opt.x))

def dice_example_p():
    """
    Example with a biased dice with a previously observed mean
    p is optimized directly, rather than the lagrange multiplier l
    """
    def entropy(p, *args):
        """
        Entropy of p
        """
        return sum(p*np.log(p))

    def sum_constraint(p, *args):
        """
        Constraint that p sum to 1
        """
        return sum(p) - 1

    def mean_constraint(p, *args):
        """
        Constraint that the weighted mean is fixed
        """
        mean = args[0]
        x = np.arange(1,7)
        return p.dot(x) - mean

    # Observed mean
    mean = 3
    # Starting guess
    p = np.ones(6)/6
    # Set constraints
    constraints = [{'type':'eq',
                    'fun':sum_constraint},
                    {'type':'eq',
                    'fun':mean_constraint,
                    'args': (mean,)}
                    ]

    opt = scipy.optimize.minimize(entropy, p, bounds=[[1e-9,None] for _ in range(6)],
            constraints=constraints)
    print(opt.x)