gamma-poissonモデルとdynamic-gamma-poissonモデル

2x2_problem.py

import math

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import gamma, poisson

# parameters
E_p0 = 0.088
g = 500
a = E_p0 * g

N_0 = 2000
N_1 = 40000

q_0 = 0.1
q_1 = 0.1

prior_dist = gamma(a, scale=1./g)

p_0 = prior_dist.rvs()
print(p_0)

xs = np.arange(0, 1.0, 0.01)
ys = []
gains = []

for x in xs:
    click_dist = poisson(p_0 * x * N_0)
    c = click_dist.rvs()

    a = a + c
    g = g + N_0 * x
    gain = N_0 * x * (E_p0 - q_0) + N_1 * max(a/g - q_1, 0)

    E_clicks = gain + q_0 * N_0 + q_1 * N_1

    ys.append(E_clicks)
    gains.append(gain)


plt.title("p_0 = {}".format(math.floor(p_0 * 1000) / 1000) + ' ' +
          "q_0 = q_1 = {}".format(q_0) + ' '
          "E[p0] = {}".format(E_p0)
          )

plt.xlabel('x')
plt.ylabel('E[#clicks]')
plt.plot(xs, ys, 'r')

plt.savefig("2x2problem.png")

dynamic_gamma_poisson.py

import math

import matplotlib.animation as animation
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import gamma, poisson

# input
p1 = 0.6
p2 = 0.2
N = 10
w = 0.95

# parameters
a = 0
g = 0

step = 50
xs = np.arange(0.01, 1.0, 0.01)

fig = plt.figure()


def update(frame):
    global a, g

    plt.cla()

    if frame <= step / 2:
        p = p1
    else:
        p = p2

    click_dist = poisson(p * N)

    c = click_dist.rvs()

    a = w * a + c
    g = w * g + N

    posterior_dist = gamma(a, scale=1./g)

    y = [posterior_dist.pdf(x) for x in xs]

    plt.plot(xs, y, "r")
    E = math.floor(a / g * 100) / 100
    plt.title('frame={}'.format(frame) + ', ' +
              'CTR={}'.format(p) + ': ' +
              'E[p_t]={}'.format(str(E).ljust(4, '0')),
              )
    plt.xlim(0, 1)
    plt.ylim(0, 15)


ani = animation.FuncAnimation(fig, update, frames=step)
ani.save('dgp.gif')

gamma_poisson.py

import math

import matplotlib.animation as animation
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import gamma, poisson

# input
p1 = 0.6
p2 = 0.2
N = 10

# parameters
a = 0
g = 0

step = 50
xs = np.arange(0.01, 1.0, 0.01)

fig = plt.figure()


def update(frame):
    global a, g

    plt.cla()

    if frame <= step / 2:
        p = p1
    else:
        p = p2

    click_dist = poisson(p * N)

    c = click_dist.rvs()

    a = a + c
    g = g + N

    posterior_dist = gamma(a, scale=1./g)

    y = [posterior_dist.pdf(x) for x in xs]

    plt.plot(xs, y, "r")
    E = math.floor(a / g * 100) / 100
    plt.title('frame={}'.format(frame) + ', ' +
              'CTR={}'.format(p) + ': ' +
              'E[p_t]={}'.format(str(E).ljust(4, '0')),
              )
    plt.xlim(0, 1)
    plt.ylim(0, 15)


ani = animation.FuncAnimation(fig, update, frames=step)
ani.save('gp.gif')

main.py

import math

import matplotlib.animation as animation
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import gamma, poisson

# input
p1 = 0.6
p2 = 0.2
N = 10
w = 0.95

# parameters
a = 0
g = 0
d_a = 0
d_g = 0

step = 50
xs = np.arange(0.01, 1.0, 0.01)

fig, ax = plt.subplots()


def update(frame):
    global a, g
    global d_a, d_g

    ax.cla()

    if frame <= step / 2:
        p = p1
    else:
        p = p2

    click_dist = poisson(p * N)

    c = click_dist.rvs()

    a = a + c
    g = g + N

    posterior_dist = gamma(a, scale=1./g)

    y = [posterior_dist.pdf(x) for x in xs]
    E = math.floor(a / g * 100) / 100

    d_a = w * d_a + c
    d_g = w * d_g + N

    d_posterior_dist = gamma(d_a, scale=1./d_g)

    d_y = [d_posterior_dist.pdf(x) for x in xs]
    d_E = math.floor(d_a / d_g * 100) / 100

    ax.set_title('frame={}'.format(frame) + ', ' +
                 'CTR={}'.format(p) + ': ' +
                 'E[p(gp)]={}'.format(str(E).ljust(4, '0')) + ' ' +
                 'E[p(dgp)]={}'.format(str(d_E).ljust(4, '0')),
                 )

    ax.set_xlim(0, 1)
    ax.set_ylim(0, 15)

    ax.plot(xs, y, color='red', label='gp')
    ax.plot(xs, d_y, color='blue', label='dpg')
    ax.legend(loc=0)


ani = animation.FuncAnimation(fig, update, frames=step)
ani.save('image.gif')