# 計算を素早くしてくれる

import numpy as np

# フレーム（エクセル）で行列の計算や集計など

import pandas as pd

# 機械学習

import seaborn as sns

# MCMC

import pystan

# グラフ

import matplotlib.pyplot as plt

# Pythonの中で図表示

%matplotlib inline

# 表示桁数の指定

%precision 3

# 複雑な統計計算

from scipy.stats import mstats

data = pd.read.csv ("ファイル名.csv", sep=";", na_values=".") # csvファイルをdataという名前をつけて読み込む。第２引数以降はオプション

data.shape　# データの数

pd.DataFrame.describe ( )

data.head ( ) # データの最初の５行目までを表示（きちんと読み込まれているか確認）

マックの場合で文字化けする場合

data = pd.read_csv("/content/drive/MyDrive/Python/フォルダ名/ファイル名.csv", encoding='utf-8')

x = data ["age"] # ｘにdataの中の"age"の列を指定

y = data ["P1"] # y にdataの中の"P1"の列を指定

np.corrcoef(x, y) # Numpyで相関係数を確認

散布図

plt.scatter(data ["P1"], ["age"]) #P1 と age の散布図

# 単回帰分析 age でP2を回帰

stan_model = """

data {

int N;

real X [N];

real Y [N];

}

parameters {

real a;

real b;

real , low = 0 sigma;

}

model {

for (n in 1 1:N) {

Y [N] ~ normal (a * X[N] + b, sigma);

}

}

"""

sm = pystan.StanModel (model_code=stan_model)

stan_data = { "N" : df.shape[0], "X" : df["age"], "Y" : df["P2"] }

fit = sm.sampling(data = stan_data, iter=2000, warmup=500, chains = 3) # 繰り返し数2,000回、ウォームアップ500回、３セット繰り返し

fit

※　確認する結果

１．Rhat値　1.0がベスト

各値を確認

fig = fit.plot()　# 事後分布とMCMCの結果

結果の平均値から予測直線を引く

a = 0.16

b =-0.89

x = np.arange(65, 100, 1)

y = a * x + b

plt.plot(x, y)

plt.scatter(df["age"], df["P2"])

stan_model = """

data {

int N ;

real X [N];

real Y [N];

int N_s;

real X_s [N_s];

}

parameters {

real a;

real b;

real < lower = 0 > sigma;

}

model {

for (n in 1:N) {

Y [n] ~ normal (a * X[n] + b, sigma);

}

}

generated quantities {

real Y_s [N_s];

for (n in 1 : N_s) {

Y_s [n] = normal_rng (a * X_s [n] + b, sigma) ;

}

}

"""

sm = pystan.StanModel(model_code=stan_model)

X_s = np.arange( 65, 100, 1)

N_s = X_s.shape[0]

stan_data = { "N" : df.shape[0], "X" : df["P2"], "Y" : df["age"] , "N_s" : N_s, "X_s" : X_s}

fit = sm.sampling(data = stan_data, iter = 2000, warmup = 500, chains = 3)

fit.extract("a") # a についてMCMCで取得したデータを抽出

ms_a = fit.extract("a")["a"]

plt.hist(ms_a)

ms_b = fit.extract("b")["b"]

df_b = pd.DataFrame([])

for i in range(65, 100,1):

df_b[i] = ms_a * i + ms_b

df_b

low_y50, high_y50 = mstats.mquantiles(df_b, [0.25, 0.75], axis = 0)

low_y95, high_y95 = mstats.mquantiles(df_b, [0.025, 0.975], axis = 0)

plt.scatter(df["age"], df["P2"])

plt.fill_between(X_s, low_y50, high_y50, alpha = 0.6, color = "darkgray")

plt.fill_between(X_s, low_y95, high_y95, alpha = 0.3, color = "gray")

a = 0.16

b =-0.89

y = a * X_s + b

plt.plot(X_s, y, color = "black")

Y_p = fit.extract("Y_s")["Y_s"]

low_y, high_y = msats.mquantiles((Y_p, [0.025,0.975], asis = 0)

plt.scatter(df["age"], df["P2"])

plt.fill_between(X_s, low_y, high_y, alpha = 0.3, color = "gray")

a = 0.16

b =-0.89

y = a * X_s + b

plt.plot(X_s, y, color = "black")

重回帰分析

必要な変数のみのデータフレームを作成

dfm = df[["P1", "age", "economy", "sex"]]

変数間の関係を図で確認

g = sns.PairGrid (dfm)

g= g.map_lower (sns.kdeplot)

g = g.map_diag (sns.distplot, kde = False)

g = g.map_upper (plt.scatter)

stan_model = """

data {

int N ;

real age [N];

real economy [N];

real sex [N];

real Y [N];

}

parameters {

real a ;

real e;

real s;

real b;

real < lower = 0 > sigma;

}

model {

real mu;

for (n in 1:N) {

mu = a * age[n] + e * economy [n] + s * sex [n] + b;

Y [n] ~ normal (mu, sigma);

}

}

"""

sm = pystan.StanModel(model_code=stan_model)

stan_data = { "N" : df.shape[0], "age" : df ["age"], "economy" : df ["economy"], "sex" : df["sex"], "Y" : df["P1"] }

fit = sm.sampling(data = stan_data, iter = 2000, warmup = 500, chains = 4)

fit

fig = fit.plot()