The probability of density function

General Usage
1

I want to test an MCMC code for a Gaussian distribution witha standard deviation of 0.1and for random values of the mean

using Distributions
using LinearAlgebra
using Plots

sigma=0.1
log_target(mu)= log((1/sigma*sqrt(pi*2^2))*exp((-(x-mu)^2)/(2*sigma^2)))
proposal(mu)=Normal(mu,sigma)
initial=0
Iteration=10000
Chain=fill(0.0,Iteration)
Chain[1]=initial
select=[]
for i in 2:Iteration
    current= Chain[i-1]
    proposed= rand(proposal(current))
            
    
    C= min(1, log_traget(proposed)/log_traget(current))

    if rand() < C
        Chain[i]= proposed
        push!(select, Chain[i])
        
        
    else
        Chain[i]=current
        
    end    
end
select

the code gives me the selected values but i know it’s not enough i want to plot the density and i have no clue how to do it i tried but it didn’t work, i’m still a beginner at mcmc

2

You plot a histogram: histogram(select)

1 Like
3

If I may, I have a simple package to plot densities and other simple fits (EasyFit):

julia> using Plots, EasyFit

julia> x = randn(1000);

julia> fit = fitdensity(x)
 
 ------------------- Density -------------

  d contains the probability of finding data points within x ± 0.0294

 ----------------------------------------- 


julia> plot(fit.x,fit.d)

julia> savefig("./plot.png")

Result:

plot

You can control some options, like the bin width:

julia> fit = fitdensity(x,step=0.5)
 
 ------------------- Density -------------

  d contains the probability of finding data points within x ± 0.25

 ----------------------------------------- 


julia> plot(fit.x,fit.d)

Result:

plot

2 Likes
4

what if i want without a histogram only the plot

5

I’m not sure I understand you. But if I do, then you may want to look into something like EasyFit or other ways of kernel density estimation. But usually in MCMC, one just plots the histograms, sometimes 2D ones to see the correlations between varibles. Have a look at GitHub - JuliaPlots/StatsPlots.jl: Statistical plotting recipes for Plots.jl

1 Like
6

Check out the MCMCChains.jl package. You can do something like

# visualize results
chn = Chains(chain, ["μ₁","μ₂","σ₁","σ₂","p"])
display(chn)
plot(chn)

and see summary information for the chain and marginal density plots.

2 Likes