How get the values of pdf function Beta

Specific Domains Statistics
1 
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 {isplaystyle athrm {B} (lpha ,eta )={rac {amma (lpha )amma (eta )}{amma (lpha +eta )}}} ?

2

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
3

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
4

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

5

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")
6

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

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

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.

8

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
9

Various beta functions are available in SpecialFunctions

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

10

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

Immagine 2023-02-21 154724
Immagine 2023-02-21 1547241500×387 33.7 KB

11

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
12

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
13

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


using QuadGK
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