Python-сообщество

# -*- coding: utf-8 -*-
from IPython.core.pylabtools import figsize
import numpy as np
import pymc as pm
from matplotlib import pyplot as plt
from matplotlib import rc
font = {'family': 'Verdana', 'weight': 'normal'}
rc('font', **font)
figsize(10, 5)
data = np.genfromtxt("v15.csv", skip_header = 1, dtype = int, delimiter=";")
data = data[np.argsort(data[:, 0])]
data = [d[1:] for d in data]
number_of_players = len(data)
number_of_checkpoints = len(data[0])
difficulty_checkpoint = pm.Uniform("dc", lower = 0, upper = 1, size = number_of_checkpoints)
skill_player = pm.Uniform("sp", lower = 0, upper = 1, size = number_of_players)
out = []
for i in range(number_of_players):
    out.append([0] * number_of_checkpoints)
@pm.deterministic
def theta(dc = difficulty_checkpoint, sp = skill_player):
#    out = np.zeros(number_of_checkpoints)
#    out = [d*s for s in sp for d in dc]
    for i in range(number_of_players):
        for j in range(number_of_checkpoints):
            out[i][j] = sp[i] * dc[j]
    return out
checkpoint1 = 0
checkpoint2 = 19
@pm.deterministic
def delta(dc = difficulty_checkpoint):
    return dc[checkpoint1] - dc[checkpoint2]
    
observation_player = np.arange(number_of_players * number_of_checkpoints)\
                        .reshape(number_of_players, number_of_checkpoints)
for i in range(number_of_players):
    for j in range(number_of_checkpoints):
        observation_player[i][j] = pm.Bernoulli("obs",
                                                p = theta[i][j],
                                                value = data[i][j],
                                                observed = True)
model = pm.Model([difficulty_checkpoint, skill_player, observation_player, delta])
mcmc = pm.MCMC(model)
mcmc.sample(5000, 100)
dc_samples = mcmc.trace("dc")[:]
sp_samples = mcmc.trace("sp")[:]
delta_samples = mcmc.trace("delta")[:]
print("\nВероятность, что КП_1 ЛЕГЧЕ КП_2: %.3f" % (delta_samples < 0).mean())
print("Вероятность, что КП_1 СЛОЖНЕЕ КП_2: %.3f" % (delta_samples > 0).mean())
plt.figure()
figsize(10, 6)
ax = plt.subplot(311)
ax.set_autoscaley_on(False)
for i in range(number_of_checkpoints):
    plt.hist(dc_samples.T[i], histtype='stepfilled', bins=30, alpha=0.33, normed=True)
    plt.legend(loc="upper right")
    plt.title(r"""Апостериорные распределения параметров""")
    plt.xlim([0, 1])
    plt.ylim([0, 10])
    plt.xlabel("Что это?")
plt.figure()
figsize(10, 6)
ax = plt.subplot(312)
ax.set_autoscaley_on(False)
for i in range(number_of_players):
    plt.hist(sp_samples.T[i], histtype='stepfilled', bins=30, alpha=0.33, normed=True)
    plt.legend(loc="upper right")
    plt.title(r"""Апостериорные распределения параметров""")
    plt.xlim([0, 1])
    plt.ylim([0, 10])
    plt.xlabel("А Что это?")
plt.figure()
figsize(10, 6)
ax = plt.subplot(313)
ax.set_autoscaley_on(False)
plt.hist(delta_samples, histtype='stepfilled', bins=30, alpha=0.85,
         label="Апостериорное распределение delta", color="#7A00A6", normed=True)
plt.ylim(0, 10)