How to set up a custom likelihood function for two variables

Questions
4

Your set up is not quite right in both cases. Put it this way, your logp should be evaluted based on the input of both the observed data and the parameters. So if you think of your logp as a function, it should has (f_1, f_2, theta_1, theta_2) as inputs. The differences between Potential and Densitydist is that, in Potential you dont distinguish between data and parameters (i.e., everything is just input to the logp), whereas in Densitydist you set up the parameterization, and evaluate on the data (the logp conditioned on the parameters).

(1) Using Potential:

def likelihood(theta_1, theta_2, f_1, f_2):
    a = m_1*m_2
    b = m_1*k_2*theta_2+m_2*k_2*theta_2+m_2*k_1*theta_1
    c = theta_1*theta_2*k_1*k_2

    # Controlled parameter
    i0 = 9

    f_1_square = (b-tt.sqrt(b**2-4*a*c))/(2*a)/(4*np.pi**2)
    f_2_square = (b+tt.sqrt(b**2-4*a*c))/(2*a)/(4*np.pi**2)

    return tt.exp(-2**(i0-2)*((f_1_square/f_1**2-1)**2+
                        (f_2_square/f_2**2-1)**2))

with pm.Model() as model:
    
    theta_1 = pm.Lognormal('theta_1', mu=mu_1, sd=sigma_1)
    theta_2 = pm.Lognormal('theta_2', mu=mu_2, sd=sigma_2) 

    like = pm.Potential('like', likelihood(theta_1, theta_2, fm_1, fm_2))
    
    trace = pm.sample()

(1) Using Densitydist:

def likelihood(theta_1, theta_2):
    
    def logp_(value):
        a = m_1*m_2
        b = m_1*k_2*theta_2+m_2*k_2*theta_2+m_2*k_1*theta_1
        c = theta_1*theta_2*k_1*k_2

        # Controlled parameter
        i0 = 9

        f_1_square = (b-tt.sqrt(b**2-4*a*c))/(2*a)/(4*np.pi**2)
        f_2_square = (b+tt.sqrt(b**2-4*a*c))/(2*a)/(4*np.pi**2)

        return tt.exp(-2**(i0-2)*((f_1_square/value[0]**2-1)**2+
                            (f_2_square/value[1]**2-1)**2))
    
    return logp_
    

with pm.Model() as model:
    
    theta_1 = pm.Lognormal('theta_1', mu=mu_1, sd=sigma_1)
    theta_2 = pm.Lognormal('theta_2', mu=mu_2, sd=sigma_2) 

    like = pm.DensityDist('like', likelihood(theta_1, theta_2), observed=[fm_1, fm_2])
    
    trace = pm.sample()

Also notices that I am sampling using the default setting. It is much better than Metropolis with find_MAP, which we now strongly discouraged.

You can see more information from other post on the discourse: Difference between DensityDist and Potential?, Problem with pymc3.DensityDist() #2846, Likelihood evaluation and DensityDist

2 Likes