「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。データ分析・語学に力点を置いています。 →現在、コンサルタントの雛になるべく、少しずつ勉強中です(※2024年1月21日改訂)。

MENU

Juliaを使ってみる(07/22):確率試行のシミュレーション

1. Juliaでの確率分布の扱い

 

1.1 単回帰モデル

 回帰係数を


\begin{aligned}
(\alpha,\beta)\sim \mathcal{N}_2(\boldsymbol{\mu},\Sigma)
\end{aligned}

と仮定し、この分布からサンプリングした上で単回帰モデル


\begin{aligned}
Y=\alpha+\beta X+\varepsilon,\ \varepsilon\sim N(0,\sigma^2)
\end{aligned}

に従う変数を作成、図示してみる。

using Distributions
using PyPlot

# 母数を生成するための分布
mu = [0.0, 0.0]
Sigma = [0.2 0.0; 0.0 0.5]

# 出力値に追加するノイズ
num_sigma = 0.6

# 入力値は事前に与える
X = [-10, -9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10]

# サンプリング回数
num_samples = 5

# 母数の標本
W = rand(MvNormal(mu, Sigma), num_samples)

# 出力値を保存するためのリスト
Ys = []

fig, axes = subplots(1,2, figsize = (8,3))
xs = range(-12,12,length = 100)
for n in 1:num_samples
    w1, w2 = W[:,n]
    
    # 母数をプロット
    axes[1].scatter(w1, w2)
    
    # 生成された関数をプロット
    f(x) = w1 * x + w2
    axes[2].plot(xs, f.(xs))
    
    # 関数からの出力値も生成する
    Y = rand.(Normal.(f.(X),num_sigma))
    push!(Ys, Y)
end

axes[1].set_xlabel("w1"), axes[1].set_ylabel("w2")
axes[1].set_title("parameters")

axes[2].set_xlabel("x"), axes[2].set_ylabel("y")
axes[2].set_title("functions")

tight_layout()

# 各標本セットをプロット
fig,axes = subplots(1, num_samples, sharey = true, figsize = (12,3))
xs = range(-12,12, length = 1000)
for n in 1:num_samples
    w1, w2 = W[:,n]
    Y = Ys[n]
    f(x) = w1 * x + w2
    axes[n].plot(xs, f.(xs))
    axes[n].scatter(X,Y)
    
    axes[n].set_xlabel("x"), axes[1].set_ylabel("y")
    axes[n].set_title("w1=$(round(w1, digits = 3)), w2 = $(round(w2, digits = 3))")
end

tight_layout() 




1.2 ロジスティック回帰モデル

 回帰係数を


\begin{aligned}
(\alpha,\beta)\sim \mathcal{N}_2(\boldsymbol{\mu},\Sigma)
\end{aligned}

と仮定し、この分布からサンプリングした上でロジスティック回帰モデル


\begin{aligned}
Z&\sim B(p),\\
p&=\exp\left(-\displaystyle{\frac{1}{1+Y}}\right),\\
Y&=\alpha+\beta X
\end{aligned}

に従う変数を作成、図示してみる。

###########################
### ロジスティック回帰モデル ###
###########################
using Distributions
using PyPlot

# シグモイド関数を定義
sig(x) = 1/(1 + exp(-x))

# 母数を生成するための分布
mu = [0.0, 0.0]
Sigma = [0.01 0.0; 0.0 0.01]

# 入力値は予め与えておく
X = [-10, -9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10]

# サンプリング回数
num_samples = 5

# 母数の標本
W = rand(MvNormal(mu, Sigma), num_samples)

# 出力値を保存するためのリスト
Ys = []

fig, axes = subplots(1, 2, figsize = (8, 3))
xs = range(-12,12, length = 100)
for n in 1:num_samples
    w1, w2 = W[:, n]
    
    # 母数をプロット
    axes[1].scatter(w1, w2)
    
    # 生成された関数をプロット
    f(x) = sig(w1 * x + w2)
    axes[2].plot(xs, f.(xs))
    
    Y = rand.(Bernoulli.(f.(X)))
    push!(Ys, Y)
end

axes[1].set_xlabel("w1"), axes[1].set_ylabel("w2")
axes[1].set_title("parameters")

axes[2].set_xlabel("x"), axes[2].set_ylabel("y")
axes[2].set_title("functions")

tight_layout()

# 各標本セットをプロット
fig,axes = subplots(1, num_samples, sharey = true, figsize = (12,3))
xs = range(-12,12, length = 1000)
for n in 1:num_samples
    w1, w2 = W[:,n]
    Y = Ys[n]
    f(x) = sig(w1 * x + w2)
    axes[n].plot(xs, f.(xs))
    axes[n].scatter(X,Y)
    
    axes[n].set_xlabel("x"), axes[1].set_ylabel("y")
    axes[n].set_title("w1=$(round(w1, digits = 3)), w2 = $(round(w2, digits = 3))")
end

tight_layout() 




1.3 ポアソン回帰モデル

 回帰係数を


\begin{aligned}
(\alpha,\beta)\sim \mathcal{N}_2(\boldsymbol{\mu},\Sigma)
\end{aligned}

と仮定し、この分布からサンプリングした上でポアソン回帰モデル


\begin{aligned}
Z&\sim Po(\lambda),\\
Z&=\exp\left(-\displaystyle{\frac{1}{1+Y}}\right),\\
Y&=\alpha+\beta X
\end{aligned}

に従う変数を作成、図示してみる。

###################
### Poisson回帰 ###
###################

using Distributions
using PyPlot

# 母数の従う分布の母数
mu = [0.0, 0.0]
Sigma = [0.01 0.0; 0.0 0.01]

# 入力値
X = [-10, -9, -8, -7, -6, -5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10]

# サンプリング回数
num_samples = 5

# 母数の標本
W = rand(MvNormal(mu, Sigma), num_samples)

Ys = []

fig,axes = subplots(1, 2, figsize = (8, 3))
xs = range(-12,12,length = 1000)
for n in 1:num_samples
    w1,w2 = W[:,n]
    
    # 母数をプロット
    axes[1].scatter(w1, w2)
    
    # 生成された関数をプロット
    f(x) = exp(w1 * x + w2)
    axes[2].plot(xs, f.(xs))
    
    # 関数からの出力値も生成
    Y = rand.(Poisson.(f.(X)))
    push!(Ys, Y)
end

axes[1].set_xlabel("w1"), axes[1].set_ylabel("w2")
axes[1].set_title("parameters")

axes[2].set_xlabel("x"), axes[2].set_ylabel("y")
axes[2].set_title("functions")

tight_layout()


fig,axes = subplots(1, num_samples, sharey = true, figsize = (12,3))
xs = range(-12,12, length = 1000)
for n in 1:num_samples
    w1, w2 = W[:,n]
    Y = Ys[n]
    f(x) = exp(w1 * x + w2)
    axes[n].plot(xs, f.(xs))
    axes[n].scatter(X,Y)
    
    axes[n].set_xlabel("x"), axes[1].set_ylabel("y")
    axes[n].set_title("w1=$(round(w1, digits = 3)), w2 = $(round(w2, digits = 3))")
end

tight_layout() 




プライバシーポリシー お問い合わせ