using Distributions

using Plots

using Random

using LaTeXStrings

Random.seed!(12345)

begin

    # Initialize

    nBurnin = 50 # burn in count

    nChains = 5  # Markov Chains

    # Define Transition Operator

    P(x,chains) = rand(Normal(0.5*x,1),1,chains) # 1x5

    nTransitions = 1000

    x = zeros(nTransitions,nChains)

    x[1,:] = randn(1,nChains) # 초기값 평균 0, 분산 1 인 정규분포에서 샘플링

end

begin

    for iT in 2:nTransitions

        x[iT,:] = P(x[iT-1],nChains)

    end

    x[end]

end

begin

    minn = min(x...)

    maxx = max(x...)

    #https://discourse.julialang.org/t/plots-jl-layout-macro-usage/11024

    l = @layout[[a ; b] c ]

    plot([nBurnin,nBurnin],[minn,maxx],linecolor=:black,linewidth=3,

            label="~Burnin",legend=:bottomright)

    p1=plot!(x[1:100,:],label="",title="First 100 Samples",titlefont=font(10))  

    plot([nBurnin,nBurnin],[minn,maxx],linecolor=:black,linewidth=3,

            label="~Burnin",legend=:bottomright)    

    p2=plot!(x,label="",title="Entire Chain",titlefont=font(10))    

    # Display Samples from stationary distribution

    samples = x[nBurnin+1:end,:]

    μ = trunc(mean(samples),digits=6)

    σ = trunc(std(samples),digits=6)

    p3=histogram(samples[:],bins=100,label="Markov Chian\nSamples",

        title=L"\mu=%$μ, \sigma=%$σ)",titlefont=font(10),alpha=0.5)

    plot(p1,p2,p3,layout=l)

end