<jupyter_start><jupyter_text>どのぐらいで収束するかを見る。どれも1万回とかで十分収束しているのでは？<jupyter_code>from scipy.stats import invgamma
from scipy.stats import norm
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
% matplotlib inline

pine = pd.read_table("pine.txt", delim_whitespace = True)
p = pine.values
ave_x = np.average(p[:, 1]).item()

mu_alpha = 3000
sigma_alpha = 10**6
mu_beta = 185
sigma_beta = 10**4
a = 3
b = 1/(2*300**2)

sample_iter = 10000
burn_in = 30
N = np.shape(p)[0]
var = np.var(p[:, 1])

def sum1(beta):
    elements = np.ones(N)
    for i in range(N):
        elements[i] = p[i, 0] - beta * (p[i, 1] - ave_x)
    return np.sum(elements)

def sum2(alpha):
    elements = np.ones(N)
    for i in range(N):
        elements[i] = (p[i, 1] - ave_x)*(p[i, 0] - alpha)
    return np.sum(elements)

def sum3(alpha, beta):
    elements = np.ones(N)
    for i in range(N):
        elements[i] = (p[i, 0] - alpha - beta*(p[i, 1] - ave_x))**2
    return np.sum(elements)

def loglike(alpha, beta, sigma):
    return N*np.log(1/(sigma*np.sqrt(2*np.pi))) - sum3(alpha, beta) / (2*sigma)

n = 20
c = 5

alphas = [3000]
betas = [185]
sigmas = [90000]

t = (1/10)**2

alphas = np.ones(sample_iter)
betas = np.ones(sample_iter)
sigmas = np.ones(sample_iter)

alphas[0] = 3000
betas[0] = 185
sigmas[0] = 90000

t = (1/10)**2

for j in range(sample_iter -1):

    location_alpha = (sigma_alpha*t*sum1(betas[-1]) + sigmas[-1]*mu_alpha) / (sigma_alpha * N*t + sigmas[-1])
    scale_alpha = np.sqrt((sigma_alpha * sigmas[-1]) / (sigma_alpha * N*t + sigmas[-1]))
    r = norm.rvs(loc = location_alpha, scale = scale_alpha)
    alphas[j+1] = r.item()

    location_beta = (sigma_beta * t * sum2(alphas[-1]) + sigmas[-1] * mu_beta) / (sigma_beta *t* var*N + sigmas[-1])
    scale_beta = np.sqrt((sigmas[-1] * sigma_beta) / (sigma_beta *t* var*N + sigmas[-1]))
    q = norm.rvs(loc = location_beta, scale = scale_beta)
    betas[j+1] = q.item()

    shape = N*t/2 + a
    invrate = 2*b / (b*t*sum3(alphas[-1], betas[-1]) + 2)
    s = invgamma.rvs(a = shape, scale = 1/ invrate)
    sigmas[j+1] = s.item()
plt.plot(betas)
plt.plot(alphas)
plt.plot(sigmas)<jupyter_output><empty_output><jupyter_text>並列化でも結果を次のルーぷに使うから無理やん<jupyter_code>from scipy.stats import invgamma
from scipy.stats import norm
from multiprocessing import Pool
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
% matplotlib inline

pine = pd.read_table("pine.txt", delim_whitespace = True)
p = pine.values
ave_x = np.average(p[:, 1]).item()

mu_alpha = 3000
sigma_alpha = 10**6
mu_beta = 185
sigma_beta = 10**4
a = 3
b = 1/(2*300**2)

# ここをいじってね。
sample_iter = 1000
burn_in = 30
m = 10
core = 3

N = np.shape(p)[0]
var = np.var(p[:, 1])

def sum1(beta):
    elements = np.ones(N)
    for i in range(N):
        elements[i] = p[i, 0] - beta * (p[i, 1] - ave_x)
    return np.sum(elements)

def sum2(alpha):
    elements = np.ones(N)
    for i in range(N):
        elements[i] = (p[i, 1] - ave_x)*(p[i, 0] - alpha)
    return np.sum(elements)

def sum3(alpha, beta):
    elements = np.ones(N)
    for i in range(N):
        elements[i] = (p[i, 0] - alpha - beta*(p[i, 1] - ave_x))**2
    return np.sum(elements)

def loglike(alpha, beta, sigma):
    return N*np.log(1/(sigma*np.sqrt(2*np.pi))) - sum3(alpha, beta) / (2*sigma)

def function(w):
    if w < 100:
        
        n = 20
        c = 3
        estimates = np.ones(n+1)
        ts = np.ones(n+1)

        for i in range(n+1):
            t = (i/n)**c
            ts[i] = t

            if i == 0:
                alphas = [3000]
                betas = [185]
                sigmas = [90000]

            else:
                alphas = [np.mean(alpha_sample).item()]
                betas = [np.mean(beta_sample).item()]
                sigmas = [np.mean(sigma_sample).item()]

            for j in range(sample_iter):

                location_alpha = (sigma_alpha*t*sum1(betas[-1]) + sigmas[-1]*mu_alpha) / (sigma_alpha * N*t + sigmas[-1])
                scale_alpha = np.sqrt((sigma_alpha * sigmas[-1]) / (sigma_alpha * N*t + sigmas[-1]))
                r = norm.rvs(loc = location_alpha, scale = scale_alpha)
                alphas.append(r.item())

                location_beta = (sigma_beta * t * sum2(alphas[-1]) + sigmas[-1] * mu_beta) / (sigma_beta * var*N + sigmas[-1])
                scale_beta = (sigmas[-1] * sigma_beta) / (sigma_beta * var*N + sigmas[-1])
                q = norm.rvs(loc = location_beta, scale = scale_beta)
                betas.append(q.item())

                shape = N*t/2 + a
                invrate = 2*b / (b*t*sum3(alphas[-1], betas[-1]) + 2)
                s = invgamma.rvs(a = shape, scale = 1/ invrate)
                sigmas.append(s.item())

            alpha_sample = alphas[burn_in:]
            beta_sample = betas[burn_in:len(betas)]
            sigma_sample = sigmas[burn_in:len(sigmas)]

            box = np.ones(len(alpha_sample))
            for k, l in enumerate(alpha_sample):
                box[k] = loglike(l, beta_sample[k], sigma_sample[k])

            estimates[i] = np.average(box)
        
        return estimates

if __name__ == '__main__':
    p = Pool(core) 
    result = p.map(function, range(m))
result<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1625428245146555, 'fraction_numerical': 0.03673391701560716, 'mean_word_length': 2.5244802146210596, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '403967', 'license': 'MIT', 'n_tokens_mistral': 2237, 'n_tokens_neox': 2074, 'n_words': 450, 'path': 'midterm.ipynb', 'repo_name': 'keiikegami/midtermtest', 'size': '5254'}