Specifying complex truncated beta distribution

235 Views Asked by At

I'd like to specify a truncated beta distribution such that the support is [0.1, 0.4], and allow for the probability density at the lower bound to be higher than very near 0.

Here's what I have so far. I realize that the parameters I specify here may not give me the distribution exactly like I described, but I'm just trying to get a truncated beta distribution to work:

import numpy as np
import pandas as pd
import plotly.express as px
from scipy import stats

def get_ab(mean, stdev):
    kappa = (mean*(1-mean)/stdev**2 - 1)
    return mean*kappa, (1-mean)*kappa

class truncated_beta(stats.rv_continuous):
    def _pdf(self, x, alpha, beta, a, b):
        return stats.beta.pdf(x, alpha, beta) / (stats.beta.cdf(b, alpha, beta) - stats.beta.cdf(a, alpha, beta))
    def _cdf(self, x, alpha, beta, a, b):
        return (stats.beta.cdf(x, alpha, beta) - stats.beta.cdf(a, alpha, beta)) / (stats.beta.cdf(b, alpha, beta) - stats.beta.cdf(a, alpha, beta))
   
def plot_beta_distr(mu, stdev, lb, ub):
    alpha_, beta_ = get_ab((mu-lb)/ub, stdev/ub)

    dist = truncated_beta(a=lb, b=ub, alpha = alpha_, beta = beta_, name='truncated_beta')

    x = np.linspace(lb, ub, 100000)
    y = dist.pdf(x, alpha = alpha_, beta = beta_, a = lb, b = ub)

    fig = px.line(x = x, y = y)
    return fig

mu = 0.02
stdev = 0.005
lb = 0.1
ub = 0.4

plot_beta_distr(mu, stdev, lb, ub)

I get an error at the truncated_beta step:

TypeError: __init__() got an unexpected keyword argument 'alpha'

Update: The reason why the above doesn't work is very likely that I'm specifying a mu that's outside of my bounds.

Here is a simpler formulation that I generated from the advice on this post.

mean = .35
std = .1

lb = 0.1
ub = 0.5

def get_ab(mu, stdev):
    kappa = (mu*(1-mu)/stdev**2 - 1)
    return mu*kappa, (1-mu)*kappa

alpha_, beta_ = get_ab((mean - lb) / ub, std/ub)

norm = stats.beta.cdf(ub, alpha_, beta_) - stats.beta.cdf(lb, alpha_, beta_)

yr=stats.uniform.rvs(size=1000000)*norm+stats.beta.cdf(lb, alpha_, beta_)
xr=stats.beta.ppf(yr, alpha_, beta_)

xr.min()  # 0.100
xr.max()  # 0.4999
xr.std()  # 0.104
xr.mean() # 0.341

Everything lines up here except the mean is off. It's clear that I'm misspecifying something.

1

There are 1 best solutions below

1
On

You are not supposed to pass alpha and beta to the truncated_beta class, since it subclasses stats.rv_continuous, and the __init__ doesn't accept those input arguments.

Rather the alpha and beta should just be passed to the call of the pdf after the truncated_beta has been initiated.

Further, the mu you specify is out of bounds as you mention, so if you try and change it to mu=0.2, it works. See code example below

import numpy as np
import plotly.express as px
from scipy import stats


def get_ab(mean, stdev):
    kappa = (mean * (1 - mean) / stdev ** 2 - 1)
    return mean * kappa, (1 - mean) * kappa


class truncated_beta(stats.rv_continuous):
    def _pdf(self, x, alpha, beta, a, b):
        return stats.beta.pdf(x, alpha, beta) / (stats.beta.cdf(b, alpha, beta) - stats.beta.cdf(a, alpha, beta))

    def _cdf(self, x, alpha, beta, a, b):
        return (stats.beta.cdf(x, alpha, beta) - stats.beta.cdf(a, alpha, beta)) / (
                    stats.beta.cdf(b, alpha, beta) - stats.beta.cdf(a, alpha, beta))


def plot_beta_distr(mu, stdev, lb, ub):
    alpha_, beta_ = get_ab((mu - lb) / ub, stdev / ub)

    dist = truncated_beta(a=lb, b=ub, name='truncated_beta')

    x = np.linspace(lb, ub, 100000)
    y = dist.pdf(x, alpha=alpha_, beta=beta_, a=lb, b=ub)

    fig = px.line(x=x, y=y)
    return fig


mu = 0.2
stdev = 0.005
lb = 0.1
ub = 0.4

my_fig = plot_beta_distr(mu, stdev, lb, ub)
my_fig.show()

Example image