Likelihood analysis is typically limited to normally distributed noise due to the difficulty of determining the probability density function of complex, high-dimensional, non-Gaussian, and anisotropic noise. This work presents Score-based LIkelihood Characterization (SLIC), a framework that resolves this issue by building a data-driven noise model using a set of noise realizations from observations. We show that the approach produces unbiased and precise likelihoods even in the presence of highly non-Gaussian correlated and spatially varying noise. We use diffusion generative models to estimate the gradient of the probability density of noise with respect to data elements. In combination with the Jacobian of the physical model of the signal, we use Langevin sampling to produce independent samples from the unbiased likelihood. We demonstrate the effectiveness of the method using real data from the Hubble Space Telescope and James Webb Space Telescope.