<jupyter_start><jupyter_text>Generate Data  Blue line is true distribution, points are data.<jupyter_code>mu0 = np.array([0.8, 0.1])
t0 = np.array([0.4, 0.1])
p0 = np.array([0.6, 0.4])

a = 1/(t0-t0*mu0)-mu0
b = 1/(t0*mu0)-1+mu0

Nbin = 100
Nsamp = 10

Nsubsamp = np.random.multinomial(Nsamp, p0)
datavals = np.array([], dtype=int)
for idx, Nss in enumerate(Nsubsamp):
    pvals = np.random.beta(a[idx],b[idx],Nss)
    datavals = np.concatenate([datavals, np.random.binomial(Nbin,pvals)])
    
data1 = dict(Nsamp=Nsamp, Nbin=Nbin, Nc=15, x=datavals, mu=np.dot(mu0,p0))

def weighted_beta(weights, mu, t):
    x = np.linspace(0,1,500)
    y = np.zeros(len(x))
    for idx in range(len(weights)):
        a = 1/(t[idx]-t[idx]*mu[idx])-mu[idx]
        b = 1/(t[idx]*mu[idx])-1+mu[idx]
        y += weights[idx] * st.beta.pdf(x, a, b)
    return x, y

x,y = weighted_beta(p0, mu0, t0)
plt.plot(x,y)
plt.plot(datavals / float(Nbin), np.zeros(Nsamp),'*');<jupyter_output><empty_output><jupyter_text>PyMC3 Model<jupyter_code>def stick_breaking(beta):
    portion_remaining = tt.concatenate([[1], tt.extra_ops.cumprod(1 - beta)[:-1]])
    return beta * portion_remaining
def scale_nu(weights, locations, mu0):
    n = weights.shape[0]
    exp_locations = tt.exp(-locations)
    
    h = pm.math.logit(mu0) - tt.dot(weights, locations)
    for i in range(2):
        tmp = pm.math.invlogit(h + locations)
        h = h - (tt.dot(weights, tmp) - mu0) / (tt.dot(weights, tmp * (1-tmp)))
        # print(h.tag.test_value)
    return pm.math.invlogit(h + locations)

with pm.Model() as model:
    alpha = pm.Gamma('alpha', 1, 1)
    q = pm.Beta('q', 1, alpha, shape=data1['Nc'])
    w = pm.Deterministic('w', stick_breaking(q))
    
    nu_star = pm.Normal('nu_star', 0, 1.9, shape=data1['Nc'], testval=np.ones((data1['Nc'])))
    t = pm.Beta('t', 1, 1, shape=data1['Nc'])
    nu = pm.Deterministic('nu', scale_nu(w, nu_star, data1['mu']))
    a_ = 1/(t * (1-nu)) - nu
    # print(a_.tag.test_value)
    obs = pm.Mixture('obs', w,
                    pm.BetaBinomial.dist(
                                        alpha = 1/(t * (1-nu)) - nu,
                                        beta = 1/(t * nu) - (1-nu),
                                        n = data1['Nbin']
                                        ),
                    observed=data1['x']
                    )
with model:
    trace = pm.sample(3000, njobs=2, tune=1000)

pm.traceplot(trace);
randlist = np.random.choice(1000, 100, replace=False)
burnin=2000
for i in np.sort(randlist):
    point = trace[burnin+i]
    x,y = weighted_beta(weights=point['w'], mu=point['nu'], t=point['t'])
    plt.plot(x, y, color='grey', alpha=.1)
    
x,y = weighted_beta(p0, mu0, t0)
plt.plot(x,y)
plt.plot(datavals / float(Nbin), np.zeros(Nsamp), '*')
plt.ylim(-.2,7);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18772053634438957, 'fraction_numerical': 0.03669724770642202, 'mean_word_length': 2.69140625, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '642873', 'license': 'MIT', 'n_tokens_mistral': 1202, 'n_tokens_neox': 1137, 'n_words': 207, 'path': 'PyMC3QnA/github_pymc3_issue_2131.ipynb', 'repo_name': 'junpenglao/Planet_Sakaar_Data_Science', 'size': '2834'}