statistics - Estimate the parameters of a random variable which is the sum of Uniform Random Variables using PyMC -

statistics - Estimate the parameters of a random variable which is the sum of Uniform Random Variables using PyMC -

i trying utilize pymc estimate parameters of model. unable understand how 1 estimates parameters of model not standard distribution perhaps sum or function of other distributions.

example: lets have info "data" generated process sum of 2 random variables x , y both drawn uniform distributions parameters (a, b) , (c,d) respectively. model using pymc , estimate parameters a, b,c, , d. able setup priors parameters not sure how specify observed variable , bind observed data.

if distribution of observed variable standard (say o) do:

obs = pm.disto(params, observed= true, value=data)

but not case. can model scenario in pymc @ ?

python code using below:

import numpy np import pymc pm # generate synthetic info = 2.0 b = 8.0 c = 6.0 d = 10.0 d1 = np.random.uniform(a, b, 100) d2 = np.random.uniform(c, d, 100) info = d1 + d2 # lets seek recover parameters. #setup priors # info observed. lets recover params p_a = pm.normal("pa", 0.0, 10.0) p_b = pm.normal("pb", 0.0, 10.0) p_c = pm.normal("pc", 0.0, 10.0) p_d = pm.normal("pd", 0.0, 10.0) p_d1 = pm.uniform("pd1", p_a, p_b) p_d2 = pm.uniform("pd2", p_c, p_d) # here confused ? # p_data = p_d1 + p_d2 # how specify p_data's value observed (the observations in "data") #todo: utilize mcmc sample , obtain traces

you can model scenario in pymc2, , in sense easy so. hard do, demonstrate solution special case of model $b-a = d-c$.

i easy because pymc2 can utilize arbitrary function info log-likelihood, using observed decorator. example:

@pm.observed def pdata(value=data, pa=pa, pb=pb, pc=pc, pd=pd): logpr = 0. # code calculate logpr info , parameters homecoming logpr

it not easy because have come code calculate logpr, , find complicated , error-prone cases sum of 2 uniforms. code in inner loop of mcmc, efficiency important.

if info single uniform unknown support, utilize decorator approach follows:

@pm.observed def pdata(value=data, pa=pa, pb=pb, pc=pc, pd=pd): homecoming pm.uniform_like(value, pa, pb)

in special case of model uniforms have same width, info likelihood proportional triangular distribution, , can write out moderate amount of pain:

@pm.observed def pdata(value=data, pa=pa, pb=pb, pc=pc, pd=pd): pd = pc + (pb - pa) # don't utilize pd value # create sure order acceptible if pb < pa or pd < pc: homecoming -np.inf x = value pr = \ np.where(x < pa+pc, 0, np.where(x <= (pa+pb+pc+pd)/2, x - (pa+pc), np.where(x <= (pb+pd), (pb+pd) - x, 0))) \ / (.5 * ((pb+pd) - (pa+pc)) * ((pb-pa) + (pd-pc))/2) homecoming np.sum(np.log(pr))

for general sum-of-two-uniforms posted, distribution trapezoidal, , writing out more work for. approximate bayesian computation (abc) promising approach situations much work compute log-likelihood explicitly.

here a notebook collecting up.

statistics pymc