Search code examples
scipylog-likelihood

Joint Likelihood in Scipy


The code below finds the joint likelihood for n (= 3 ) variables from a normal distribution with known sigma. I'm trying to extend this to finding the joint likelihood for n variables from a location-scale-t distribution with three unknown parameters (location, scale, and df), but am struggling with how exactly to proceed here. Is there a helpful function for this? At the moment I am attempting to do this manually by specifying multiple parameter domains...

n = 3; x = stats.norm(loc=0,scale=1).rvs(n, random_state = 10); # set seed
parameter_domain = np.linspace(-10,10,1001);
likelihood = stats.norm.pdf((x[:,np.newaxis])*np.ones(parameter_domain.shape),
                            loc=parameter_domain1, scale=1)#.prod(axis=0) # sigma=1 known
n = 3
df_true = 15, loc_true = 10, scale_true = 2 
x = stats.t(df = df_true, loc = loc_true, scale = scale_true).rvs(n, random_state = 10) # Data from location-scale-t distn.
df_parameter_domain = np.linspace(1, 100, 100)#np.linspace(-10,10,1001);
loc_parameter_domain = np.linspace(-10,10,1001)
scale_parameter_domain = np.linspace(0, 10, 1001)

likelihood = stats.t.pdf((x[:,np.newaxis])*np.ones(parameter_domain.shape),
                            loc=df_parameter_domain, scale=scale_parameter_domain, df = df_parameter_domain)#.prod(axis=0) # no variables known

Solution

  • The following code contains a function to calculate the likelihood you need.

    from scipy import stats
    
    n = 3
    df_true = 15
    loc_true = 10
    scale_true = 2
    x = stats.t(df=df_true, loc=loc_true, scale=scale_true).rvs(n, random_state=10)
    
    
    def likelihood_fun(loc, scale, df):
        return stats.t.pdf(x, loc, scale, df).prod()
    
    
    likelihood = likelihood_fun(loc=loc_true, scale=scale_true, df=df_true)
    

    In this example, the likelihood is evaluated at the true parameters, but they can be replaced by any value. The likelihood has a 3D vector as input and the return value is a 1D scalar. The function likelihood_fun can be vectorized but it will be complicated and possibly not very useful. If you intend to maximize it to find a maximum likelihood estimator, this implementation will do the job. If you want to plot it, you will need to fix one or two parameters to display it respectively as a 2D surface or as a 1D function.