In general we maximize a function
$$ L(\theta; x_1, \ldots, x_n) = \prod_{i=1}^n f(x_i \mid \theta) $$
where $f$ is probability density function if the underlying distribution is continuous, and a probability mass function (with summation instead of product) if the distribution is discrete.
How do we specify the likelihood function if the underlying distribution is a mixture between a continuous and a discrete distribution, with the weights on each depending on $\theta$ ?
I admit to puzzling over this question for quite some time earlier in my career. One way I convinced myself of the answer was to take an extremely practical, applied view of the situation, a view that recognizes no measurement is perfect. Let's see where that might lead.
The point of this exercise is to expose the assumptions that might be needed to justify the somewhat glib mixing of densities and probabilities in expressions for likelihoods. I will therefore highlight such assumptions wherever they are introduced. It turns out quite a few are needed, but they're pretty mild and cover every application I have encountered (which obviously will be limited, but still includes quite a few).
The problem concerns a mixed distribution $F,$ one that is neither absolutely continuous nor singular. Lebesgue's Decomposition Theorem permits us to view such a distribution as a mixture of an absolutely continuous one (which by definition has a density function $f_a$) and a singular ("discrete") one, which has a probability mass function $f_d.$ (I'm going to ignore the possibility that a third, continuous but not absolutely continuous component, may be present. Those who use such models tend to know what they're doing and usually have all the technical skills to justify them.)
When $F = F_\theta$ is a member of a parametric family of distributions, we may write
$$F_\theta(x) = F_{a\theta}(x) + F_{d\theta}(x) = \int_{\infty}^x f_a(t;\theta)\mathrm{d}t + \sum_{t \le x} f_d(t;\theta).$$
(The sum is at most countable, of course.) Here, $f_a(\,;\theta)$ is a probability density function multiplied by some mixture coefficient $\lambda(\theta)$ and $f_d(\,;\theta)$ is a probability mass function multiplied by $1-\lambda(\theta).$
Let's interpret any observation $x_i$ in an iid dataset $X=(x_1,x_2,\ldots, x_n)$ as "really" meaning we have certain knowledge that a hypothetical true underlying value $y_i$ lies in an interval $(x_i-\delta_i, x_i+\epsilon_i]$ surrounding $x_i,$ but otherwise have no information about $y_i.$ Assuming we know all the deltas and epsilons, this no longer presents any problems for constructing a likelihood because everything can be expressed in terms of probabilities:
$$\mathcal{L}(X;\theta) = \prod_i \left(F_\theta(x_i + \epsilon_i) - F_\theta(x_i - \delta_i)\right).$$
If the support of $F_{d\theta}$ has no condensation points at any $x_i,$ its contribution to the probability will reduce to at most a single term provided the epsilons and deltas are made sufficiently small: there will be no contribution when $x_i$ is not in its support.
If we assume $f_a(\,;\theta)$ is Lipschitz continuous at all the data values, then uniformly in the sizes of the epsilons and deltas we may approximate the absolutely continuous part of $F_\theta(x_i)$ as
$$F_{a\theta}(x_i + \epsilon_i) - F_{a\theta}(x_i - \delta_i) = f_a(x_i;\theta)(\epsilon_i + \delta_i) + o(|\epsilon_i + \delta_i|).$$
The uniformity of this approximation means that as we take all the epsilons and deltas to grow small, all the $o()$ terms also grow small. Consequently there is a vanishingly small value $\epsilon(\theta)\gt 0,$ governed by the contributions of all these error terms, for which
$$\eqalign{ \mathcal{L}(X;\theta) &= \prod_i \left(f_a(x_i;\theta)(\epsilon_i + \delta_i) + o(|\epsilon_i + \delta_i|) + f_d(x_i;\theta)\right)\\ &= \prod_i \left(f_a(x_i;\theta)(\epsilon_i + \delta_i) + f_d(x_i;\theta)\right)\ + \ o(\epsilon(\theta)). }$$
This is still a little messy, but it shows where we're going. In the case of censored data, usually just one part of each term in the product will be nonzero, because these models typically assume that the support of the singular part of the distribution is disjoint from the upport of the continuous part, no matter what the parameter $\theta$ might be. (Specifically: $f_d(x) \ne 0$ implies $F_a(x+\epsilon)-F_a(x-\epsilon) = o(\epsilon).$) That permits us to break the product into two parts and we can factor the contributions from all the intervals out of the continuous part:
$$\mathcal{L}(X;\theta) = \left(\prod_{i=1}^k (\epsilon_i + \delta_i) \right)\prod_{i=1}^k f_a(x_i;\theta) \ \prod_{i=k+1}^n f_d(x_i;\theta).$$
(Without any loss of generality I have indexed the data so that $x_i, i=1, 2, \ldots, k$ contribute to the continuous part and otherwise $x_i, i=k+1, k+2, \ldots, n$ contribute to the singular part of the likelihood.)
This expression now makes it plain that
Since the interval widths $\epsilon_i+\delta_i$ are fixed, they do not contribute to the likelihood (which is defined only up to some positive constant multiple).
Accordingly, we may work with the expression
$$\mathcal{L}(X;\theta) = \prod_{i=1}^k f_a(x_i;\theta) \ \prod_{i=k+1}^n f_d(x_i;\theta)$$
when constructing likelihood ratios or maximizing the likelihood. The beauty of this result is that we never need to know the sizes of the finite intervals that are used in this derivation: the epsilons and deltas drop right out. We only need to know that we can make them small enough for the likelihood expression we actually work with to be an adequate approximation to the likelihood expression we would use if we did know the interval sizes.
This question is an extremely important foundational problem in likelihood analysis, and also a very subtle and difficult one, so I'm quite surprised at some of the superficial answers it is receiving in the comments.
In any case, in this answer I am just going to add one small point to whuber's excellent answer (which I think is the correct approach to this problem). That point is that likelihood functions in this context come from density functions over a mixed dominating measure, and this leads to the interesting property that we can scale the relative sizes of the likelihood function arbitrarily over the continuous and discrete parts and we still have a valid likelihood function. This gives rise to an obvious question of how we can implement likelihood techniques when there is no unique likelihood function.
Illustrating this point requires some preliminary presentation on the sampling density as a Radon-Nikodym derivative of the probability measure, so please bear with me. I will first show how to obtain a density function for a mixed dominating measure and then I will show why this leads to the ability to scale the continuous and discrete parts of the likelihood at will. Finally, I will discuss the implications of this issue for likelihood-based analysis and give my opinion on its resolution. I think this is essentially resolved by the method whuber presents in his answer, but it would need to be extended in the direction I have discussed in the comments to that answer, so as to ensure that each point in the support of the discrete part ignores the continuous part at that point.
Expressing the density using a dominating measure: The standard approach to dealing with mixed densities for real random variables is to use Lebesgue measure $\lambda_\text{LEB}$ as the dominating measure for the continuous part and counting measure $\lambda_\text{COUNT}$ (over some specified countable set $\mathcal{D} \subset \mathbb{R}$) as the dominating measure for the discrete part. This leads to the Radon-Nikodym derivative defined by:
$$\mathbb{P}(X \in \mathcal{A} | \theta) = \int \limits_\mathcal{A} f(x | \theta) \ d \lambda_\text{LEB}(x) + \int \limits_\mathcal{A} p(x | \theta) \ d\lambda_\text{COUNT}(x).$$
(Note that the latter integral degenerates down to a sum over the elements $x \in \mathcal{A} \cap \mathcal{D}$. We write it here as an integral to make the similarity between the two terms clearer.) One can use a single density by taking the measure $\lambda_* \equiv \lambda_\text{LEB} + \lambda_\text{COUNT}$ and setting:
$$f_*(x | \theta) \equiv \mathbb{I}(x \notin \mathcal{D}) \cdot f(x | \theta) + \mathbb{I}(x \in \mathcal{D}) \cdot p(x | \theta).$$
Using $\lambda_*$ as the dominating measure, we then have the following expression for the probability of interest:
$$\mathbb{P}(X \in \mathcal{A} | \theta) = \int \limits_\mathcal{A} f_*(x | \theta) \ d \lambda_*(x).$$
This shows that the function $f_*$ is a valid Radon-Nikodym derivative of the probability measure on $X$, so it is a valid density for this random variable. Since it depends on $x$ and $\theta$ we can then define a valid likelihood function $L_x^*(\theta) \propto f_*(x | \theta)$ by holding $x$ fixed and treating this as a function of $\theta$.
Effect of scaling the dominating measures: Now that we understand the extraction of a density from a dominating measure, this leads to a strange property where we can scale the relative sizes of the likelihood over the continuous and discrete parts and we still have a valid likelihood function. If we now use the dominating measure $\lambda_{**} \equiv \alpha \cdot \lambda_\text{LEB} + \beta \cdot \lambda_\text{COUNT}$ for some positive constants $\alpha > 0$ and $\beta > 0$ then we now get the corresponding Radon-Nikodym density:
$$f_{**}(x | \theta) \equiv \frac{\mathbb{I}(x \notin \mathcal{D})}{\alpha} \cdot f(x | \theta) + \frac{\mathbb{I}(x \in \mathcal{D})}{\beta} \cdot p(x | \theta).$$
Using $\lambda_{**}$ as the dominating measure, we then have the following expression for the probability of interest:
$$\mathbb{P}(X \in \mathcal{A} | \theta) = \int \limits_\mathcal{A} f_{**}(x | \theta) \ d \lambda_{**}(x).$$
As in the above case, we can define a valid likelihood function $L_x^{**}(\theta) \propto f_{**}(x | \theta)$ by holding $x$ fixed and treating this as a function of $\theta$. You can see that the freedom to vary $\alpha$ and $\beta$ now gives us freedom to scale the relative sizes of the continuous and discrete parts in the likelihood function as much as we want, and still have a valid likelihood function (albeit with respect to a different dominating measure, with corresponding scaling of parts).
This particular result is just part of the more general result that every likelihood function is defined with respect to some (implicit) underlying dominating measure, and there is no unique likelihood function that can be defined irrespective of this underlying measure.$^\dagger$ Nevertheless, in this particular case we see that it is still based on a dominating measure that is a combination of Lebesgue measure and counting measure, so we have not really monkeyed with the measure very much. Since there is no objective justification for forming the dominating measure from equal weightings of Lebesgue measure and counting measure, the implication of this is that there is no objective justification for the relative scaling for the continuous and discrete parts of the likelihood function.
Implications for likelihood analysis: This might seem to put us in a bit of a quandary. We can arbitrarily scale the discrete and continuous parts of the likelihood function up or down in relative size and still have just as reasonable a claim to this being a valid likelihood function. Fortunately, this problem can be solved by recognising that the scaling constants will come out of the likelihood function in the same way as illustrated in whuber's answer. That is, if we have $x_1,...,x_k \notin \mathcal{D}$ and $x_{k+1},...,x_n \in \mathcal{D}$ we will get:
$$\begin{equation} \begin{aligned} L_\mathbb{x}^{**}(\theta) = \prod_{i=1}^n L_{x_i}^{**}(\theta) &= \prod_{i=1}^n f_{**}(x_i | \theta) \\[12pt] &= \Bigg( \prod_{i=1}^k \frac{1}{\alpha} \cdot f(x_i | \theta) \Bigg) \Bigg( \prod_{i=k+1}^n \frac{1}{\beta} \cdot p(x_i | \theta) \Bigg) \\[12pt] &= \frac{1}{\alpha^k \beta^{n-k}} \Bigg( \prod_{i=1}^k f(x_i | \theta) \Bigg) \Bigg( \prod_{i=k+1}^n p(x_i | \theta) \Bigg) \\[12pt] &= \frac{1}{\alpha^k \beta^{n-k}} \prod_{i=1}^n f_{*}(x_i | \theta) \\[12pt] &\propto \prod_{i=1}^n f_{*}(x_i | \theta) \\[12pt] &= \prod_{i=1}^n L_{x_i}^{*}(\theta) \\[12pt] &= L_\mathbb{x}^{*}(\theta). \\[12pt] \end{aligned} \end{equation}$$
This shows that the scaling properties of the dominating measure only affect the likelihood function through a scaling constant that can be ignored in standard MLE problems. Note that in my treatment of this problem, this useful property has occurred as a direct result of the fact that the sampling density is defined in a way that ignores the continuous density when we are in the support of the discrete part. (This differs from whuber's answer, where he allows for a combination of these parts. I think this might actually lead to some difficult problems; see my comments to that answer.)
$^\dagger$ This result is not confined to mixed cases. Even in simple cases with continuous or discrete random variables, if you vary the underlying dominating measure it will give a corresponding variation in the Radon-Nikodym derivative, which then leads to a different likelihood function.
The likelihood function $\ell(\theta|\mathbf{x})$ is the density of the data at the observed value $\mathbf{x}$ expressed as a function of $\theta$ $$\ell(\theta|\mathbf{x})=f(\mathbf{x}|\theta)$$ This density is defined for every (acceptable) value of $\theta$ almost everywhere over the support of $\mathbf{x}$, $\mathfrak{X}$, against a particular measure over $\mathfrak{X}$ that does not depend on $\theta$. For any parametric family, there should exist such a dominating measure across all $\theta$'s, hence a density, hence a likelihood.
Here is a relevant excerpt from the Wikipedia entry on likelihood functions (stress is mine):
In measure-theoretic probability theory, the density function is defined as the Radon-Nikodym derivative of the probability distribution relative to a dominating measure. This provides a likelihood function for any probability model with all distributions, whether discrete, absolutely continuous, a mixture or something else. (Likelihoods will be comparable, e.g., for parameter estimation, only if they are Radon–Nikodym derivatives with respect to the same dominating measure.)
One example where this occurs, that is, likelihood given by a probability model of mixed continuous/discrete type, is with censored data. For an example see Weighted normal errors regression with censoring .
In general this can be formulated using measure theory. Then assume a statistical model with a model function $f(x;\theta)$ which is a Radon-Nikodym derivative with respect to a common measure $\lambda$ (which should not depend on the parameter $\theta$). Then the likelihood function based on an independent sample $x_1, x_2, \dotsc, x_n$ is $\prod_i f(x_i;\theta)$. This is really the same in continuous, discrete and mixed cases.
One simple example could be modeling of daily rainfall. That could be zero, with positive probability, or positive. So for the dominating measure $\lambda$ we could use the sum of Lebesgue measure on $(0,\infty)$ and an atom at zero.
