# How get the values of pdf function Beta

``````using Distributions
using StatsPlots
plot(Beta(1, 1), ylim=(0, 5), size=(400, 300), label=false, xlabel="x", ylabel="P(x)", title="Beta distribution shapes")

``````

How does plot get the values of the Beta pdfunction to plot the graph.
In other words, if I wanted to calculate the value of f(x; a,b) for x=0.73 how could I do it, apart from going to the definition I find on Wikipedia.

And for the calculation of Beta(a,b), how do you do it?

Perhaps in this way ?

If you want the PDF, just ask for it : `pdf(Beta(2,3),0.73)` should do the trick.

More generally, refer to the documentation of Distributions.jl : type `?pdf` in your repl

1 Like

`StatsPlots` is basically doing something like

``````julia> using Plots, Distributions

julia> plot(x -> pdf(Beta(1, 1), x), minimum(support(Beta(1,1))):0.01:maximum(support(Beta(1,1))))
``````

when you plot a distribution (while probably being a bit smarter than this about automatic x and y lims looking at the plot that the above produces…)

1 Like

Tanks!
And for the evaluation of Beta(5,3)?

Tanks!!

I tried like this and the graph is identical

``````plot(x -> pdf(Beta(10, 3), x), minimum(support(Beta(1,1))):0.01:maximum(support(Beta(1,1))), ylim=(0, 5), size=(400, 300), label=false, xlabel="x", ylabel="P(x)", title="Beta distribution shapes")
``````

Actually `support` doesn’t even need the specific beta distribution:

``````julia> support(Beta)
RealInterval{Float64}(0.0, 1.0)
``````

What do you mean by the evaluation of `Beta(5,3)` ? Do you want the beta function ? in which case,

``````using SpecialFunctions
?beta
``````

should guide you.

I’m sorry I wrote wrong. B(a,b) not Beta(a,b)
That’s what I meant …

``````a,b= 5, 3
B(a,b)=integral[0,1] u^(a-1)*(1-u)^(b-1) du
``````

Various beta functions are available in `SpecialFunctions`

https://specialfunctions.juliamath.org/stable/functions_list/#SpecialFunctions.beta

Ok.
looks like what I assume in the previous message.
But it doesn’t look as good as in the preview… I don’t know why

I’m not sure I’m following - what is your issue? The title suggests you are interested in a pdf which suggests you are talking about a distribution, hence the initial answers you got relating to the Beta probability distribution. The beta function shares a name with and is used in the pdf of the Beta distribution but is not the same thing.

``````julia> using Distributions, SpecialFunctions

julia> pdf_myBeta(α, β, x) = 1/beta(α, β) * x^(α-1) * (1-x)^(β-1)
pdf_myBeta (generic function with 1 method)

julia> pdf_myBeta(5, 3, 0.5)
1.640624999999998

julia> pdf(Beta(5, 3), 0.5)
1.6406249999999987
``````
1 Like

Nothing more than that.
I wanted to understand how plot could get the values to plot a particular beta pdf.
And since the constant B(a,b) is used in the definition of pdfbeta, I wanted to understand where and how it is calculated.
Doing @edit Beta(5,3) takes me inside the beta.jl module where I couldn’t find how the calculations of the pdf and the constant B(a,b) are done.
Now following your(s) instructions I found that it is calculated (inside SpecialFunctions.jl module, using C libraries

``````beta(x, y)
"""
Euler integral of the first kind ``\\operatorname{B}(x,y) = \\Gamma(x)\\Gamma(y)/\\Gamma(x+y)``.
"""
function beta(a::Number, b::Number)
lab, sign = logabsbeta(a, b)
return sign*exp(lab)
end

if Base.MPFR.version() >= v"4.0.0"
function beta(y::BigFloat, x::BigFloat)
z = BigFloat()
ccall((:mpfr_beta, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Int32), z, y, x, ROUNDING_MODE[])
return z
end
end

``````

a bit for fun I tried to get the results also in the following ways

``````
a,b=5,3
integral, err = quadgk(t -> t^(a-1)*(1-t)^(b-1), 0, 1, rtol=1e-20)
# (0.009523809523809523, 0.0)
using SpecialFunctions
beta(5,3) # 0.009523809523809535

using Integrals
f(t,p) = p*t^4*(1-t)^2
p=1
prob = IntegralProblem(f, 0, 1, p)
sol = solve(prob) # u: 0.009523809523809523

setprecision(60, base=10) # use 60-digit arithmetic

integral, err = quadgk(t -> t^(a-1)*(1-t)^(b-1), big"0.0", big"1.0", rtol=1e-50)
# (0.0095238095238095238095238095238095238095238095238095238095238286, 9.7234613716580339174126000840314441259222690136268236454704947e-63)

beta(5,3) ≈ integral #true
beta(5,3) ≈  factorial(4)*factorial(2)/factorial(7) #true
beta(5,3) ≈ gamma(5)*gamma(3)/gamma(8) #true

# "≈" can be typed by \approx<tab>
``````
1 Like