Pluto notebook: Animation of Gaussian process regression

GP.jl

### A Pluto.jl notebook ###
# v0.14.2

using Markdown
using InteractiveUtils

# ╔═╡ 9dcadaad-5e5e-494f-931d-1ab9269ed8c1
using Plots, LinearAlgebra

# ╔═╡ 077d65e8-a167-4527-9bb8-8fa9eb614c46
let
	n = 50; X = range(0, 6π; length=n); Y = sin.(X) + 0.2randn(n)
	L = 600; x = range(-0.5π, 7π; length=L)
	g(x₁, x₂; a=1.0, b=2.0, c=0.2) = a*exp(-(x₁ - x₂)^2/(2*b^2)) + c*(x₁ == x₂)
	anim = @animate for k in [fill(1, 5); 1:n; fill(n, 10)]
		U, V = @view(X[1:k]), @view(Y[1:k])
		A = g.(x, U')/g.(U, U'); μ = A*V; S = g.(x, x') - A*g.(U, x'); d = .√diag(S)
		scatter(U, V; label="", msw=0, xlim=extrema(x), ylim=(-3, 3))
		plot!(x, [μ μ]; label="", fillrange=[μ-2d μ+2d], color=:red, fillalpha=0.1)
	end; gif(anim, "GP.gif"; fps=10)
end

# ╔═╡ 469ff1fd-b8c8-4678-8d68-f79aa4f0554d
let
	# サンプル生成
	n = 50
	X = range(0, 6π; length=n)
	Y = sin.(X) + 0.2randn(n)
	
	# 予測区間
	L = 600
	x = range(-0.5π, 7π; length=L)
	
	# モデル
	g(x₁, x₂; a=1.0, b=2.0, c=0.2) = a*exp(-(x₁ - x₂)^2/(2*b^2)) + c*(x₁ == x₂)
	
	# ガウス過程回帰
	anim = @animate for k in [fill(1, 5); 1:n; fill(n, 10)]
		# サンプル生成
		U, V = @view(X[1:k]), @view(Y[1:k])
		
		# ガウス過程回帰 (予測分布は平均μ, 分散共分散行列Sの多変量正規分布になる)
		KXX, KXx, KxX, Kxx = (g.(u, v') for u in (U, x) for v in (U, x))
		A = KxX/KXX
		μ = A*V
		S = Kxx - A*KXx
		d = .√diag(S)
		
		# プロット
		scatter(U, V; label="", msw=0, xlim=extrema(x), ylim=(-2.5, 2.5))
		plot!(x, [μ μ]; label="", fillrange=[μ-2d μ+2d], color=:red, fillalpha=0.1)
	end
	gif(anim, "GP.gif"; fps=10)
end

# ╔═╡ Cell order:
# ╠═9dcadaad-5e5e-494f-931d-1ab9269ed8c1
# ╠═077d65e8-a167-4527-9bb8-8fa9eb614c46
# ╠═469ff1fd-b8c8-4678-8d68-f79aa4f0554d