今回やりたいこと
確率試行のシミュレーションをJuliaで扱ってみる。
を参照する。
1. Juliaでの確率分布の扱い
1.1 単回帰モデル
回帰係数を
と仮定し、この分布からサンプリングした上で単回帰モデル
に従う変数を作成、図示してみる。
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 ロジスティック回帰モデル
回帰係数を
と仮定し、この分布からサンプリングした上でロジスティック回帰モデル
に従う変数を作成、図示してみる。
########################### ### ロジスティック回帰モデル ### ########################### 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 ポアソン回帰モデル
回帰係数を
と仮定し、この分布からサンプリングした上でポアソン回帰モデル
に従う変数を作成、図示してみる。
################### ### 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()
