Thursday, 26 January 2023

Specifying complex truncated beta distribution

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.



from Specifying complex truncated beta distribution

No comments:

Post a Comment