# The probability of density function

I want to test an MCMC code for a Gaussian distribution witha standard deviation of `0.1`and 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

You plot a histogram: `histogram(select)`

1 Like

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:

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:

2 Likes

what if i want without a histogram only the plot

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

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