Kinetic modeling of alcohol consumption

TL;DR

This paper models alcohol consumption distributions using kinetic theory, explaining their shape through microscopic agent interactions and deriving a generalized Gamma distribution as the equilibrium.

Contribution

It introduces a kinetic model based on prospect theory to explain the emergence of Gamma and Weibull distributions in alcohol consumption data.

Findings

Alcohol consumption distributions can be modeled by a generalized Gamma distribution.

Microscopic agent interactions determine the macroscopic distribution shape.

The model clarifies the meaning of parameters in the equilibrium distribution.

Abstract

In most countries, alcohol consumption distributions have been shown to possess universal features. Their unimodal right-skewed shape is usually modeled in terms of the Lognormal distribution, which is easy to fit, test, and modify. However, empirical distributions often deviate considerably from the Lognormal model, and both Gamma and Weibull distributions appear to better describe the survey data. In this paper we explain the appearance of these distributions by means of classical methods of kinetic theory of multi-agent systems. The microscopic variation of alcohol consumption of agents around a universal \emph{social} accepted value of consumption, is built up introducing as main criterion for consumption a suitable value function in the spirit of the prospect theory of Kahneman and Twersky. The mathematical properties of the value function then determine the unique macroscopic equilibrium which results to be a generalized Gamma distribution. The modeling of the microscopic kinetic interaction allows to clarify the meaning of the various parameters characterizing the generalized Gamma equilibrium.

Tables

Table 1. Table 4.1. Expectation and variance of the Fokker-Planck and of the Boltzmann model for different values of the scaling factor ϵ italic-ϵ \epsilon and the shape factor δ 𝛿 \delta for the Test 1.

$\gamma =1$, ${\overline{w}}_L=6$		$E\left[g_{\infty }^B\left(w\right)\right]$	$Var\left(g_{\infty }^B\left(w\right)\right)$
$\delta =0.1$	$ϵ=0.9$	2.66	24.05
	$ϵ=0.5$	3.03	18.51
	$ϵ=0.1$	3.32	13.56
	$ϵ=0.001$	3.36	18.21
	FP	3.37	18.22
$\delta =0.3$	$ϵ=0.9$	2.19	19.58
	$ϵ=0.5$	2.52	14.93
	$ϵ=0.1$	2.81	14.93
	$ϵ=0.001$	2.85	15.05
	FP	2.85	14.93
$\delta =0.5$	$ϵ=0.9$	1.60	14.45
	$ϵ=0.5$	1.93	11.31
	$ϵ=0.1$	2.20	11.52
	$ϵ=0.001$	2.23	11.61
	FP	2.25	11.65
$\delta =0.8$	$ϵ=0.9$	0.38	3.57
	$ϵ=0.5$	0.77	4.70
	$ϵ=0.1$	1.02	5.64
	$ϵ=0.001$	1.06	5.80
	FP	1.14	6.17

Table 2. Table 4.2. Expectation and variance of the Fokker-Planck and of the Boltzmann model for different values of the scaling factor ϵ italic-ϵ \epsilon and the shape factor δ 𝛿 \delta for the Test 2.

$\gamma =2$, ${\overline{w}}_L=3$		$E\left[g_{\infty }^B\left(w\right)\right]$	$Var\left(g_{\infty }^B\left(w\right)\right)$
$\delta =0.1$	$ϵ=0.9$	2.27	6.92
	$ϵ=0.5$	2.13	3.77
	$ϵ=0.1$	2.24	3.37
	$ϵ=0.001$	2.25	3.27
	FP	2.25	3.09
$\delta =0.3$	$ϵ=0.9$	2.10	6.09
	$ϵ=0.5$	1.98	3.37
	$ϵ=0.1$	2.09	3.05
	$ϵ=0.001$	2.11	3.01
	FP	2.11	2.89
$\delta =0.5$	$ϵ=0.9$	1.91	5.34
	$ϵ=0.5$	1.82	3.04
	$ϵ=0.1$	1.94	2.78
	$ϵ=0.001$	1.96	2.79
	FP	1.96	2.70
$\delta =0.8$	$ϵ=0.9$	1.60	4.33
	$ϵ=0.5$	1.54	2.59
	$ϵ=0.1$	1.67	2.44
	$ϵ=0.001$	1.69	2.41
	FP	1.70	2.43

Equations

$$\frac{\partial g(w,t)}{\partial t}=\frac{\lambda }{2}\frac{{\partial }^2}{\partial w^2}\left(w^{2}g(w,t)\right)+\frac{\mu }{2}\frac{\partial }{\partial w}\left[\frac{1}{\delta }\left({\left(\frac{w}{{\bar{w}}_L}\right)}^{\delta }-1\right)wg(w,t)\right].$$

$$g_{\infty }(w)=g_{\infty }({\bar{w}}_L){\left(\frac{w}{{\bar{w}}_L}\right)}^{\gamma /\delta -2}\exp \left\{-\frac{\gamma }{{\delta }^2}\left({\left(\frac{w}{{\bar{w}}_L}\right)}^{\delta }-1\right)\right\},$$

$$g_{\infty }(w)=g_{\infty }({\bar{w}}_L){\left(\frac{w}{{\bar{w}}_L}\right)}^{-2}\exp \left\{-\frac{\gamma }{2}\left({\left(\log \frac{w}{{\bar{w}}_L}\right)}^2\right)\right\},$$

$${\int }_{D}f(w,t)dw$$

$${\int }_{{\mathbb{R}}_{+}}f(w,t)dw=1.$$

$$w_{\ast }=w-\Phi (w/{\bar{w}}_L)w+\eta w.$$

$${\Phi }^{\prime \prime }(s)\le 00\le s<1,{\Phi }^{\prime \prime }(s)\ge 0s>1.$$

$$-\Phi \left(1-s\right)>\Phi \left(1+s\right),$$

$${\Phi }^{\prime }\left(1+s\right)<{\Phi }^{\prime }\left(1-s\right).$$

$${\Phi }_0^{ϵ}(s)=\mu \frac{s^{ϵ}-1}{s^{ϵ}+1},s\ge 0.$$

$$-\mu \le {\Phi }_0^{ϵ}(s)\le \mu .$$

$${\Phi }_1^{ϵ}(s)=\mu \frac{e^{ϵ(s-1)}-1}{e^{ϵ(s-1)}+1},s\ge 0,$$

$${\Phi }_1^{ϵ}(s)\ge -\mu \frac{1-e^{-ϵ}}{1+e^{-ϵ}}.$$

$${\Phi }_{\delta }^{ϵ}(s)=\mu \frac{e^{ϵ(s^{\delta }-1)/\delta }-1}{e^{ϵ(s^{\delta }-1)/\delta }+1},s\ge 0,$$

$${\Phi }_{\delta }^{\overline{ϵ}}(0)=-D(\delta )=-\mu \frac{1-e^{-\overline{ϵ}/\delta }}{1+e^{-\overline{ϵ}/\delta }},$$

$$\frac{d}{dt}\int_{\mathbb{R}_{+}}\varphi(w)\,f(w,t)\,dx=\frac{1}{\tau}\Big{\langle}\int_{\mathbb{R}_{+}}\bigl{(}\varphi(w_{*})-\varphi(w)\bigr{)}f(w,t)\,dw\Big{\rangle}.$$

$$\eta \to \sqrt{ϵ}\eta .$$

$$m(t)={\int }_{{\mathbb{R}}_{+}}wf(w,t)dw$$

$$\frac{d}{dt}m(t)=ϵ{\bar{w}}_L\frac{1}{\tau }{\int }_{{\mathbb{R}}_{+}}\frac{1}{ϵ}{\Phi }_{\delta }^{ϵ}\left(w/{\bar{w}}_L\right)\frac{w}{{\bar{w}}_L}f(w,t)dw.$$

$$\frac{1}{ϵ}{\Phi }_{\delta }^{ϵ}\left(s\right)\le \mu (s-1),$$

$$\frac{1}{ϵ}{\Phi }_{\delta }^{ϵ}\left(s\right)s\ge -\mu .$$

$$\tau \to ϵ\tau ,$$

$$\frac{d}{dt}{\int }_{{\mathbb{R}}_{+}}w{f}_{ϵ}(w,t)dw={\bar{w}}_L\frac{1}{\tau }{\int }_{{\mathbb{R}}_{+}}\frac{1}{ϵ}{\Phi }_{\delta }^{ϵ}\left(w/{\bar{w}}_L\right)\frac{w}{{\bar{w}}_L}{f}_{ϵ}(w,t)dw,$$

$$A_{\delta ,ϵ}\left(w\right)=\frac{1}{ϵ}{\Phi }_{\delta }^{ϵ}\left(w/{\bar{w}}_L\right)\to \frac{\mu }{2\delta }\left({\left(\frac{w}{{\bar{w}}_L}\right)}^{\delta }-1\right).$$

$$⟨w_{\ast }-w⟩=-ϵA_{\delta ,ϵ}(w)w;⟨(w_{\ast }-w{)}^2⟩=\left({ϵ}^2{A}_{\delta ,ϵ}^2(w)+ϵ\lambda \right)w^2.$$

$$⟨\phi (w_{\ast })-\phi (w)⟩=ϵ\left(-{\phi }^{\prime }(w)w\frac{\mu }{2\delta }\left({\left(\frac{w}{{\bar{w}}_L}\right)}^{\delta }-1\right)+\frac{\lambda }{2}{\phi }^{\prime \prime }(w)w^2\right)+R_{ϵ}(w),$$

$$\frac{1}{ϵ}{R}_{ϵ}(w)\to 0$$

$$\frac{d}{dt}{\int }_{{\mathbb{R}}_{+}}\phi (w)f_{ϵ}(w,t)dw=$$

$${\int }_{{\mathbb{R}}_{+}}\left(-{\phi }^{\prime }(w)w\frac{\mu }{2\delta }\left({\left(\frac{w}{{\bar{w}}_L}\right)}^{\delta }-1\right)+\frac{\lambda }{2}{\phi }^{\prime \prime }(w)w^2\right)f_{ϵ}(w,t)dw+\frac{1}{ϵ}{\mathcal{R}}_{ϵ}(t),$$

$${\mathcal{R}}_{ϵ}(t)={\int }_{{\mathbb{R}}_{+}}{R}_{ϵ}(w)f_{ϵ}(w,t)dw.$$

$$\frac{d}{dt}{\int }_{{\mathbb{R}}_{+}}\phi (w)g(w,t)dw=$$

Full text

Department of Mathematics and Informatics of the University of Ferrara e.mail: giacomo.dimarco@unife.it

Department of Mathematics of the University of Pavia, and IMATI CNR, Italy. e.mail: giuseppe.toscani@unipv.it.

Keywords. Alcohol consumption; Empirical distribution; Generalized Gamma distribution, Log-Normal distribution, Weibull distribution, Kinetic models; Fokker–Planck equations.

1. Introduction

Alcohol consumption is a component cause for a number of diseases, which makes studies on it very relevant from both medical and psychological point of views. Information about this social phenomenon requires to know at best its statistical distribution, to analyze the effect of the resulting alcohol consumption distribution on the estimation of the alcohol Population-Attributable Fractions, and to characterize the chosen alcohol consumption distribution by exploring if there is a global relationship within the distributions [43, 59].

Despite the increasing number of papers devoted to a precise fitting of the statistical distribution of the phenomenon of alcohol consumption and to the social motivations about it, the precise connections leading from the human behavior of alcohol consumers to the observed statistical distributions have been left largely unexplored, and a satisfactory mathematical modeling is still missing.

We started to care of this problem after the reading of two interesting research papers recently published into the journal Population Health Metrics [43, 59], whose conclusions were subsequently considered as a starting point for further statistical studies (cf. [57] and the references therein). In these papers the analysis of the fitting of real data about alcohol consumption in a huge number of countries, led the authors to the conclusion that, among various probability distributions often used in this context, Gamma and Weibull distributions appeared to furnish a better fitting with respect to the Lognormal distribution, usually used in this context, starting from the work of Ledermann [44], who proposed it as a reasonable model for the consumption problem.

Concerning the mathematical modeling of alcohol consumption by resorting to the methods of statistical physics, the results of the analysis presented in [43, 59] introduce into the problem of characterizing its distribution a number of challenging questions. Given a system of agents characterized in time by their level of alcohol consumption, and given a certain observed statistical distribution of their consumption, can we describe at the level of the single agents the stylized facts which lead at a macroscopic level to this distribution? And moreover, can we characterize at the level of single agents the differences which produce a variety of steady probability distributions for the same problem?

As far as the Lognormal distribution is concerned, the kinetic modeling of phenomena concerned with interacting multi-agent systems and leading to this distribution has been recently proposed in [39], starting from the results of [38].

Among probability distributions on $\mathbb{R}_{+}$, Lognormal distribution has a leading role [1, 26]. The list of phenomena which fit Lognormal distribution in natural sciences is quite long, as documented the exhaustive review paper by Limpert, Stahel and Abbt [48]. Moreover, the success of Lognormal distribution in the study of random variations that occur in the data from these phenomena that show more or less skewed probability distributions, is favored by its features, which render it easy to fit, test and modify.

In [39], it was pointed out that a relevant number of phenomena involving measurable quantities of a population and fitting Lognormal distribution comes from social sciences and economics, areas where it can be reasonably assumed that the appearance of this distribution is a consequence of a certain human behavior.

Among others, a good fitting has been observed while looking at the distribution of body weight [18], at women’s age at first marriage [58], at drivers behavior [40], or, from the economic world, when looking at consumption in a western society [5], at the size of cities [6], and at call-center service times [17].

These phenomena were modelled in [39] by resorting to the methodology of classical kinetic theory, and the Lognormal distribution (the macroscopic equilibrium of the multi-agent system) appears as a direct consequence of the details of the microscopic interactions of agents. This idea is fully coherent with the classical kinetic theory of rarefied gases, where the formation of a Maxwellian equilibrium in the spatially uniform Boltzmann equation is closely related to the conservation laws of the elastic microscopic binary collisions between molecules [19, 37].

The analysis of [38, 39] suggests a way to modify the microscopic interaction of the variation of the alcohol consumption of agents, that leads to the class of distributions fitted in [43, 59]. As in [38, 39], the elementary single change will be expressed in terms of a suitable value function, that meets the same properties as the one leading to the Lognormal distribution, built up in the spirit of the prospect theory of Kahneman and Twersky [41, 42]. As usual in the kinetic setting [30], the microscopic interaction generates a linear kinetic equation for the density of agents, that will be subsequently studied in the asymptotic regime of grazing interactions [66]. The grazing regime describes the situation in which a single interaction produces only a very small change of the consumption variable. In this regime, the variation of density of the alcohol consumption of the agent’s system is driven by a partial differential equation of Fokker–Planck type [30].

The precise mathematical model is the following. Let us denote by $g=g(w,t)$ the density of agents which have a (positive) level of consumption equal to $w$ at time $t\geq 0$. Then, the density $g$ is the solution of a linear Fokker–Planck equation with variable coefficients of diffusion and drift, given by

[TABLE]

In (1.1) $\bar{w}_{L}$ represents the critical (maximal) level of individual alcohol consumption in the chosen time unity, while $1-\delta$, with $0\leq\delta\leq 1$ is linked to the level of psychophysical addiction induced by alcohol . Moreover, $\lambda$ and $\mu$ are positive constants closely related to the randomness and intensity of the phenomenon under study, satisfying further bounds imposed by the physical conditions of the microscopic interaction. In particular, if $\gamma=\mu/\lambda$, $\gamma>\delta$, and the equilibrium density of the Fokker–Planck equation (1.1) is given by the (non normalized) generalized Gamma density [47, 63]

[TABLE]

Moreover, for any given initial distribution $g_{0}(w)$ of agents, convergence to equilibrium can be shown to hold exponentially fast in time with explicit rate [55]. Note that the steady distribution (1.2) includes as particular cases the Gamma distribution, and in this case $\delta=1$, as well as the Weibull distribution, and in this other case $\gamma=\delta(1+\delta)$. Moreover, a (non normalized) Lognormal distribution

[TABLE]

is included into (1.2) and it is obtained from (1.2) by taking the limit $\delta\to 0$. Hence, one can conjecture that the best fitting of the alcohol consumption distribution is achieved by suitably choosing the three parameters characterizing the generalized Gamma distribution, namely $\gamma,\delta$ and $\bar{w}_{L}$.

It is remarkable that the fitting analysis of [43, 59], was limited to the choice of particular cases of the generalized Gamma distribution (1.2), obtained by fixing precise values of the parameter $\delta$. Hence, in [43, 59] the fitted distributions are characterized by the presence of only two parameters.

In the kinetic description leading to the steady state (1.2), one of these parameters is given by $\bar{w}_{L}$, representing the critical (maximal) level of individual alcohol consumption. Following Ledermann [44], this parameter can be fixed a priori, thus reducing the two-parameter distributions considered in [43, 59] to one-parameter distributions. This generates a new way of fitting that follows along the lines indicated by Ledermann. In the case of the Lognormal distribution the description of the variation of the Lognormal distribution with respect to the unique unknown parameter is given by the so-called Ledermann curve [44], which describes the percentage of consumers drinking above a fixed level. Analogous results can be derived in the other cases.

The mathematical methodology of this paper takes essential advantage from previous attempts to model social and economic phenomena in a multi-agent system. Starting from economics, where the genesis of the formation of Pareto curves in wealth distribution of western countries was studied in various aspects [20, 21, 22, 25, 27, 35, 60], this field of research moved to social sciences, where the investigation of opinion formation played a leading rule [10, 11, 12, 13, 14, 15, 16, 23, 28, 31, 32, 33, 34, 62, 64].

Recently, some of these investigations introduced into the microscopic description of the interactions specific behavioral aspects of agents. To our knowledge, this idea has been first proposed by Pareschi and Maldarella in [51]. There the authors, to investigate the price formation of a good in a multi-agent market, introduced a kinetic model for a multi-agent system consisting of two different trader populations, playing different rules of trading, and possibly changing their point of view. The kinetic description was inspired by the microscopic Lux–Marchesi model [49, 50] (cf. also [45, 46]). In [51] the trading rules of agents were assumed to depend on the opinion of traders through a kinetic model of opinion formation recently introduced in [64]. Also, psychological and behavioral components of the agents, like the way they interact with each other and perceive risks, were taken into account. This has been done by resorting, in agreement with the prospect theory by Kahneman and Twersky [41, 42], to interactions containing a suitable value function.

The analysis of [51] enlightened in a kinetic description the role of human behaviors [7, 8, 9], pioneered in [68]. These aspects have been recently considered in [38], where the introduction of a suitable value function justified at a microscopic level the mechanism of formation of the service time distribution in a call center.

The leading idea in [38] is based on a general principle which can be easily verified in a number of social activities of agents in which one identifies the possibility of a certain addiction.

The paper is organized as follows. In Section 2, we shall introduce a linear kinetic model for a multi-agent system, in which agents are characterized in terms of their level of alcohol consumption that can be measured in terms of a fixed unit of measure, and subject to microscopic interactions which describe the microscopic rate of change of the value of their consumption, according to the previous general principle. The relevant mechanism of the microscopic interaction is based on a suitable value function, in the spirit of the analysis of Kahneman and Twersky [41, 42], which reproduces at best the human behavior in this situation.

In Section 3 we will show that in a suitable asymptotic procedure (hereafter called grazing limit) the solution to the kinetic model tends towards the solution of the Fokker-Planck type equation (1.1). Similar asymptotic analysis was performed in [24, 29] for a kinetic model for the distribution of wealth in a simple market economy subject to microscopic binary trades in presence of risk, showing formation of steady states with Pareto tails, in [65] on kinetic equations for price formation, and in [64] in the context of opinion formation in presence of self-thinking. A general view about this asymptotic passage from kinetic equations based on general interactions towards Fokker–Planck type equations can be found in [30]. Other relationships of this asymptotic procedure with the classical problem of the grazing collision limit of the Boltzmann equation in kinetic theory of rarefied gases have been recently enlightened in [37].

Once the Fokker–Planck equation (1.1) has been derived, in Section (4) we will investigate at a numerical level the relationships between the solutions of the linear kinetic model of Boltzmann type and its grazing asymptotics in the form of the Fokker–Planck equation, at different values of the vanishing $\epsilon$-parameter leading to the asymptotic itself. The numerical computations will confirm the perfect agreement between the stationary solutions of the two models, even for values of the $\epsilon$-parameter of the order $10^{-1}$. This justifies, at least numerically, the role of the chosen value functions in the asymptotic description of the alcohol consumption phenomenon.

2. A linear Boltzmann equation

2.1. Modeling assumptions

Following the well-consolidated approach furnished by the kinetic theory [30, 54, 56], the statistical description of alcohol consumption in a society will be described by resorting to a linear Boltzmann-type equation in which the unknown is the probability distribution of alcohol consumption of the agent system. Since we want to fit the results obtained by this approach with the various dataset coming from real situations, the kinetic model is built up to respect some basic hypothesis we enumerate below.

First of all, given a population of agents, the population needs to be considered homogeneous with respect to alcohol consumption. In this case, to have a homogenous population, it is essential to restrict it with respect to some characteristics, like age, sex and social class [44, 61]. Once the homogeneity assumption is satisfied, agents in the system are considered indistinguishable [56]. This means that an agent’s state at any instant of time $t\geq 0$ is completely characterized by the level $w\geq 0$ of its alcohol consumption. This level can be easily and uniformly measured with respect to some unit. For example, $w$ indicates the mean value of the daily consumption of alcohol in grams of agents. The unknown is the density (or distribution function) $f=f(w,t)$, where $w\in\mathbb{R}_{+}$ and the time $t\geq 0$, and the target is to study its time evolution towards a certain equilibrium.

The precise meaning of the density $f$ is the following. Given the system of agents to study, and given an interval or a more complex sub-domain $D\subseteq\mathbb{R}_{+}$, the integral

[TABLE]

represents the number of individuals with are characterized by a level $w\in D$ of consumption at time $t>0$. It is assumed that the density function is normalized to one, that is

[TABLE]

The change in time of the density is due to the fact that agents continuously upgrade the level $w$ of their alcohol consumption in time by some action. To maintain the connection with classical kinetic theory of rarefied gases, we will always refer to a single upgrade of the level measure as a microscopic interaction.

2.2. The value function in alcohol consumption

Similarly to the problems treated in [38, 39], alcohol consumption deeply depends on both addiction and social attitudes, that can be easily translated by saying that, at least over a certain level, agents likely tend to increase the value $w$ by interactions, while manifest a certain resistance to decrease it.

For the situation under study we can identify (and fix) for the system of agents two universal values. The first value, denoted by $\bar{w}$, is the social barrier, and identifies the mean amount of alcohol consumption relative to the homogenous class, which is normally accepted by the social rules of the underlying society. The second value, denoted by $\bar{w}_{L}$, with $\bar{w}_{L}>\bar{w}$, is the threshold value that it would be better not to exceed to enter into addiction. Consequently, the human tendency to increase the value $w$ by interactions has to be coupled with the existence of this limit value $\bar{w}_{L}$ which serious health reasons suggest not to cross. Similar behaviors have been recently modeled in [38, 39] by resorting to the prospect theory introduced by Kahneman and Twersky in their pioneering paper [41], devoted to describe decision under risk.

Following [38, 39], we will describe the microscopic variation of alcohol consumption of agents in the form

[TABLE]

According to (2.4), in a single interaction the value $w$ of alcohol consumption can vary by two different motivations. The first one determines a variation of a quantity proportional to $w$ in which the coefficient $\Phi(\cdot)$, which can assume assume both positive and negative values, is a function characterizing the predictable behavior of agents. The second one, still proportional to $w$, takes into account a certain amount of human unpredictability. The random variable $\eta$ describes the unpredictable variations, which in the mean are assumed negligible, and in any case are not so significant to produce a sensible change of the value $w$. Hence, the consumption of alcohol can be both increasing and decreasing, and the mean intensity of this variation is fully determined by the function $\Phi$.

The function $\Phi$ plays the role of the value function in the prospect theory of Kahneman and Twersky [41], and contains the mathematical translation of the expected human behavior in the phenomenon under consideration.

In their pioneering theory [41, 42], Kahneman and Twersky, by resorting to various situation concerned with decision theory, identified the main properties characterizing a value function. This function has to be positive and concave above the reference value $1$ ($w>\bar{w}_{L}$), while negative and convex below ($w<\bar{w}_{L}$). Hence, in terms of the variable $s=w/\bar{w}_{L}$, provided that the value function is enough regular in $s$

[TABLE]

Thus, as far as financial transactions are considered, the value function is generally concave for gains and commonly convex for losses. Also, by noticing that the aversiveness of symmetric fair bets generally increases with the size of the stack, they assumed that the value function satisfies the further properties

[TABLE]

and

[TABLE]

Hence, properties of value functions are well defined for deviations from the reference point ($s=1$ in our case), and the value function for losses is steeper than the value function for gains.

In [38] we adapted the shape of the value function given in [41] to the problem of the characterization of the distribution of the length of a service time in a call center. A careful reading of [38] allows to recognize that the properties of the value function do not have a universal validity, since its precise form depends on the main aspects of the human behavior relative to the phenomenon under study. At difference with the choice of Kahneman and Twersky, in [38] a suitable value function was identified in the increasing concave function

[TABLE]

In (2.8) $0<\mu<1$ and $0<\epsilon<1$ are suitable constants characterizing the agents behavior. In particular, the value $\mu$ denotes the maximal amount of variation of service time that agents will be able to obtain in a single interaction. Note indeed that the value function $\Phi_{0}(s)$ is such that

[TABLE]

The function in (2.8) satisfies properties (2.6) and (2.7), but in consequence of its uniform concavity does not satisfy assumption (2.5). Note that in the service time problem studied in [38], property (2.6) characterizes in a proper way the behavior of agents which have to do a service by respecting a certain time. Indeed, given two agents starting at the same distance from the limit value $\bar{w}_{L}$ from below and above, it will be easier from the agent starting below to move closer to the optimal time, than for the agent starting above.

Going back to the form of the microscopic interaction (2.4), the minus sign in front of the value function $\Phi$ reflects the tendency of agents to increase their value $w$ of alcohol consumption when $w<\bar{w}_{L}$, while they would like to decrease consumption if $w>\bar{w}_{L}$.

As shown in [38, 39], the choice of the value function (2.8) leads, as explained in details in the forthcoming Section 3, to a stationary distribution in the form of a Lognormal distribution. Hence, this choice would be in full agreement with the data fitting of alcohol consumption proposed by Ledermann in 1956 [44], choice which is still used in present times (cf. the recent paper [53] and the references therein).

On the other hand, as briefly discussed in the introduction, the recent fitting analysis of alcohol consumption in [43, 59] seems to corroborate the hypothesis that, with respect to the fitting provided through a Lognormal distribution, a better fitting can be obtained by resorting to Gamma and Weibull distributions. These results pose a question about the correctness of the shape of the value function (2.8), which could fail to represent at best the human behavior relative to alcohol consumption.

One of the points in which the shape of the value function $\Phi_{0}^{\epsilon}(s)$, defined in (2.8), could be not so well-adapted to the description of alcohol consumers, is its behavior for small values. Indeed, the first derivative of $\Phi_{0}(s)$ in $s=0$ is blowing up. Hence, agents with a very low level of alcohol consumption are pushed very easily towards a higher consumption, action that does not agree with the common behavior of low drinkers. Indeed, agents which are characterized by a level $w\ll w_{L}$, usually change their alcohol habits very slowly. Further, the lower bound in (2.9) does not depend on the value of the critical consumption value $\bar{w}_{L}$. Hence, in terms of $\Phi_{0}(s)$ the maximal allowed increase in alcohol consumption in a single interaction does not depend on the distance of the point $s$ from the reference point, which does not appear to be supported from previous studies [61].

For these reasons, for the description of alcohol consumption, the value function

[TABLE]

seems more adapted. Similarly to (2.8), in (2.10) $0<\mu<1$ and $0<\epsilon<1$ are suitable constants characterizing the intensity of interactions of agents. At difference with the value function (2.8), the function (2.10) satisfies all properties originally proposed by Kahneman and Twersky in [41], including the property to be positive and concave above the reference value ($s>1$), while negative and convex below ($s<1$). Moreover, at difference with (2.8), while the upper bound does not depend on $\epsilon$, the lower bound depends on it, and

[TABLE]

This property implies that, for small values of the parameter $\epsilon$, the maximal amount of increase of the level of alcohol consumption is of the order of $\epsilon$.

Together with (2.10) and (2.8) we consider the class of value functions

[TABLE]

characterized by the presence of a parameter $0<\delta<1$, that interpolate between $\Phi_{0}^{\epsilon}(s)$ and $\Phi_{1}^{\epsilon}(s)$. The meaning of the $\delta$ parameter can be obtained by noticing that, for fixed values of $s$ and $\epsilon$, the value function (2.12), considered as function of $\delta$, is increasing. Hence, for values of $s<1$, which are such that $\Phi_{\delta}^{\epsilon}(s)<0$, the function $|\Phi_{\delta}^{\epsilon}(s)|$ decreases with $\delta$. For example, at the point $s=0$, given a fixed value $\bar{\epsilon}$ of $\epsilon$, we obtains

[TABLE]

where the function $D(\delta)$ is decreasing from the value $\mu$, assumed for $\delta\to 0$. Looking at the action of the value function in the elementary interaction (2.4), it appears evident that the interval in which the value function plays a major role is the interval $0\leq s\leq 1$, where we are taking into account agents that are increasing their consume of alcohol, starting from $w=0$. Hence, the value of the parameter $\delta$ is essentially linked to the greater or lesser degree of dependence exerted by the consumption of alcohol relative to agents characterized by a medium-low consumption. The previous discussion clarifies that small values of the parameter $\delta$ correspond to high values of increasing consumption and viceversa. In other words, values of $\delta$ close to zero represent the situation in which it is easier to increase the value of alcohol consumption starting from a low level.

In Figure 2.1 we report the shapes of $\Phi_{\delta}^{\epsilon}(s)$ for different choices of the parameter $\delta$.

As proven in the appendix, these functions satisfy properties (2.6) and (2.7), which render them suitable as value functions.

2.3. The kinetic model

Given the elementary interaction (2.4), and fixed a value function $\Phi$, the study of the time-evolution of the distribution of the level of alcohol consumption follows by resorting to kinetic collision-like models [19, 56]. The variation of the density $f(w,t)$ obeys to a linear Boltzmann-like equation, fruitfully written in weak form. The weak form corresponds to say that the solution $f(w,t)$ satisfies, for all smooth functions $\varphi(w)$ (the observable quantities) the integro-differential equation

[TABLE]

In (2.14) the expectation $\langle\cdot\rangle$ takes into account the presence of the random parameter $\eta$ in (2.4). The positive constant $\tau$ measures the interaction frequency.

The right-hand side of equation (2.14) represents the variation in density between agents that modify their value from $w$ to $w_{*}$ (loss term with negative sign) and agents that change their value from $w_{*}$ to $w$ (gain term with positive sign).

In reason of the nonlinearity (in the variable $w$) of the interaction (2.4), one easily verifies that equation (2.14) possesses only one conserved quantity, which is obtained by setting $\varphi=1$. This conservation law implies that the solution to (2.14) departing from a probability density function remains a probability density for all subsequent times $t>0$. Note that, provided the initial density has a certain number of moments bounded, their evolution is difficult to follow analytically. However, as shown in [38], one can obtain time-dependent upper bounds for the mean value, and for higher moments. These bounds depend of the analogous bounds for the value function, as given by (2.9). These bounds do not change by changing the value of the parameter $\delta$.

3. Grazing limit and the Fokker-Planck equation

3.1. The grazing limit

The linear kinetic equation (2.14) describes the evolution of the density consequent to interactions of type (2.4), and it is valid for any choice of the parameters $\delta,\mu$ and $\lambda$. In real situations, however, it happens that a single interaction determines only an extremely small change of the value $w$. This situation is well-known in kinetic theory of rarefied gases, where interactions of this type are called grazing collisions [56, 66]. In the value functions (2.12) the smallness assumption requires to fix $\epsilon\ll 1$. At the same time [30], the balance of this smallness with the random part is achieved by setting

[TABLE]

In this way the scaling assumptions allow to retain the effect of all parameters in (2.4) in the limit procedure. An exhaustive discussion on these scaling assumptions can be found in [30] (cf. also [39] for analogous computations in the case of the Lognormal distribution). For these reasons, we address the interested reader to these review papers for details. In what follows we only outline the differences in the asymptotic procedure induced by the new value functions (2.12). Using (3.15), we write the evolution of the mean value

[TABLE]

in the form

[TABLE]

It follows from Lagrange theorem that, if $s\geq 1$

[TABLE]

and, for $s\leq 1$

[TABLE]

These bounds are identical to the ones obtained in [38] for $\Phi_{0}^{\epsilon}$. Hence, provided that the second moment of the initial density is bounded, the scaling (3.15) is such that, for any given fixed time $t>0$, the consequent variation of the mean value $m(t)$ is small with $\epsilon$ small. To observe an evolution of the average value independent of $\epsilon$, it is enough to resort to a scaling of the frequency $\tau$. If we set

[TABLE]

and $f_{\epsilon}(w,t)$ will denote the density corresponding to the scaled interaction and frequency, then the evolution of the average value for $f_{\epsilon}(w,t)$ satisfies

[TABLE]

namely a computable evolution law for the average value of $f$, which remains bounded even in the limit $\epsilon\to 0$, since pointwise

[TABLE]

As explained in [30], the motivation leading to the scaling of the frequency of interactions is clear. Since for $\epsilon\ll 1$ the microscopic interactions produce a very small change of the value $w$, a finite variation of the mean density can be observed only if agents in the system undergo a huge number of interactions in a fixed period of time. By using the same scaling one can easily observe a well-defined evolution equation for the second moment of $f_{\epsilon}(w,t)$, which will be well-defined also in the limit $\epsilon\to 0$ (cf. the analysis in [30]).

The next step is to clarify why the linear kinetic model (2.14) modifies with $\epsilon$ to produce in the limit $\epsilon\to 0$ the Fokker–Planck type equation (1.1) [30]. Given a smooth function $\varphi(w)$, and a collision of type (2.4) that produces a small variation of the difference $w_{*}-w$, using (3.15) and (3.16) one obtains

[TABLE]

Therefore, equating powers of $\epsilon$, it holds

[TABLE]

where the remainder term $R_{\epsilon}(w)$, for a suitable $0\leq\theta\leq 1$ is such that

[TABLE]

as $\epsilon\to 0$. Therefore, the time variation of the (smooth) observable quantity $\varphi(w)$ satisfies

[TABLE]

where $\mathcal{R}_{\epsilon}(t)$ denotes the integral remainder term

[TABLE]

Letting $\epsilon\to 0$, shows that in consequence of the scalings (3.15) and (3.16) the weak form of the kinetic model (2.14) is well approximated by the weak form of a linear Fokker–Planck equation (with variable coefficients)

[TABLE]

In (LABEL:m-13) the density function $g(w,t)$ coincides with the limit, as $\epsilon\to 0$, of the density $f_{\epsilon}(w,t)$ [30]. Provided the boundary terms produced by the integration by parts vanish, equation (LABEL:m-13) coincides with the weak form of the Fokker–Planck equation

[TABLE]

Equation (3.19) describes the evolution of the distribution density $g(w,t)$ of the level of alcohol consumption $w\in\mathbb{R}_{+}$, in the limit of the grazing interactions. As often happens with Fokker-Planck type equations, the steady state density can be explicitly evaluated, and it results to be a generalized Gamma density, with parameters linked to the details of the microscopic interaction (2.4).

3.2. Steady states are generalized Gamma densities

Let us set $\gamma=\mu/\lambda$, and suppose that $\gamma>\delta$. In this case the stationary distribution of the Fokker–Planck equation (3.19) is an integrable function which can be found explicitly by solving the first order differential equation

[TABLE]

Indeed, solving (3.20) with respect to $h(w)=w^{2}g(w)$ by separation of variables allows to conclude that the unique solutions to (3.20) are the functions

[TABLE]

If we fix the mass of the steady state (3.21) equal to one, the consequent probability density is a generalized Gamma. The generalized Gamma is characterized in terms of a shape $\kappa>0$, a scale parameter $\theta>0$, and an exponent $\delta>0$, and reads [47, 63]

[TABLE]

The shape and the scale paremeter of the equilibrium state (3.21) are given by

[TABLE]

The limit $\delta\to 0$ in the Fokker–Planck equation (3.19) corresponds to the drift term induced by the value function (2.8). In this case, the equilibrium distribution (3.21) takes the form of a Lognormal density [38]

[TABLE]

where

[TABLE]

Note that for all values $\delta>0$ the moments are expressed in terms of the parameters $\bar{w}_{L}$, $\delta$, $\lambda$, and $\mu$, denoting respectively the limit level of alcohol consumption, the level of addiction the variance $\lambda$ of the random effects and the values $\delta$ and $\mu$ characterizing the value function $\phi_{\delta}^{\epsilon}$ defined in (2.12). The condition $\gamma>\delta$, that guarantees the existence of a steady state density, is achieved for small values of $\delta$, or for sufficiently large values of $\gamma$. Since $\gamma=\mu/\lambda$, the quotient between the maximal amount of level of change of alcohol consumption and the variance of the random change, the condition is satisfied when the random effects are very small.

It is interesting to remark that, in the case of the gamma distribution, the mean value of the equilibrium density (3.21) is always less than $\bar{w}_{L}$. In this case in fact

[TABLE]

Remark 3.1**.**

The analysis of the possible variety of steady states densities generated by the value functions (2.12) allows to conclude that the best fitting of the level of alcohol consumption could be obtained by a generalized Gamma distribution. Indeed, the three cases treated in [43, 59] are included as particular cases of generalized Gamma distributions. The data in [43] enlighten the fact that, in the Weibull case, the characteristic exponent $\delta\geq 1/2$. This helps to understand why the Lognormal distribution, which is obtained by pushing $\delta$ to zero, probably fails to give a good fitting of the level of alcohol consumption. It is remarkable that in [57], by noticing that the generalized Gamma nests the Gamma, Weibull and Lognormal distributions as special cases, it has been proposed as useful test distribution for the problem under investigation, since, if all other distributions are rejected, the generalized Gamma results minimally restrictive and less susceptible to misspecification bias [52].**

Remark 3.2**.**

Once the shape of the alcohol consumption distribution has been characterized, the next step would be to find the parameters of the generalized Gamma that provide the best possible fit of real datasets. However, the parameter estimation of the generalized Gamma distribution is a very difficult problem. This topic was dealt with in many papers (cf. the results in [2, 36, 52, 67] and the list of references therein) but since almost all of the proposed methods are quite complex, this topic is still an open area, that needs to be better studied.**

4. Numerical experiments

The main interest related to the Fokker–Planck equation (3.19) is the analytical description of the equilibrium, that in the present case is given by the generalized Gamma distribution (3.21). Except in some simple cases [3, 4], kinetic models of Boltzmann type do not allow similar results. On the other hand, the kinetic equation (2.14) is the right-one to study, since it contains all the microscopic details that are apparently lost in the grazing limit procedure. For this reason, we report in this Section several numerical experiments that describe convergence to the steady state (identified by the generalized Gamma distribution in the limit $\epsilon\to 0$ in the Boltzmann dynamics (2.4) for various values of the parameter $\epsilon$. In details, our study proceeds as follows. We consider two different values for the two parameters $\gamma$ and $\bar{w}_{L}$ of the model. The first parameter fixes the ratio between the intensity of the value function $\Phi_{\delta}^{\epsilon}$ and the variance of a random variation in the behavior of the agents. In the Boltzmann dynamics, we consider a uniform random variable with zero mean and variance $1/12$. The second parameter $\bar{w}_{L}$ measures the threshold value that it would be better not to exceed to enter into alcohol addiction. The values considered in the experiments are $[\gamma,\bar{w}_{L}]=[1,6]$ and $[\gamma,\bar{w}_{L}]=[2,3]$.

For these two set of parameters, we study the behavior of the microscopic dynamics detailed in (2.4) for different values of the parameter $\delta$ in the value function. This parameter is of paramount importance, since it modifies the shape of the steady states in the Fokker–Planck description. Depending on its value, from the generalized Gamma distribution we obtain a Gamma density ($\delta=1$), a Weibull density ($\delta(1+\delta)=\gamma$) or a Lognormal distribution ($\delta\rightarrow 0$) as equilibrium states.

In order to show that the Fokker-Planck dynamics can be seen as a particular case of the Boltzmann one, we place ourselves in the same scaling in which the Fokker-Planck equation (3.19) has been analytically derived. For this reason, we consider in the value function (2.12) a decreasing sequence of values for the scaling factor $\epsilon$ ranging from $0.9$ to $10^{-3}$. Correspondingly we scale the effects of the random part, according to (3.15), and the interaction frequency according to (3.16) to have an $\epsilon$-free evolution of the mean. In Figure 4.2 we show the different equilibrium states obtained with the Boltzmann dynamics for the different values of the scaling factor $\epsilon$. The Figure 4.3 reports the same results of Figure 4.2 with magnification around the origin to better appreciate the differences between the solutions. Figure 4.4 and 4.5 show the same dynamics where $\gamma$ and $\bar{w}_{L}$ have been changed to respectively $2$ and $3$, namely by diminishing the threshold and increasing the variability.

The curves have been obtained by using $10^{3}$ agents and by averaging the steady state solution over $10^{5}$ realizations. In each image it is also reported the steady state of the Fokker-Planck model (3.21) and the initial data, which are the same for both models. The initial value has been considered uniformly distributed in $[0.5,2.5]$ for all different situations analyzed. From top left to bottom right the shape factor $\delta$ takes the values $0.1$, $0.3$, $0.5$ and $0.8$ which causes the equilibrium state to completely change its shape. For all tested situations an extremely good agreement between the Boltzmann and the Fokker–Planck models is obtained from approximately $\epsilon=0.001$. However, very similar steady state profiles are found from $\epsilon=0.1$, which fully describes a kinetic evolution of Boltzmann type. This interesting behavior can be observed also from the results collected in the Tables 4.1 and 4.2, which report the expected value and the variance obtained with the kinetic Boltzmann model for the different values of the parameter $\epsilon$ and the shape factor $\delta$.

5. Appendix

In this short appendix we verify that, for any given constant $0<\delta<1$ and $\epsilon>0$, the value function $\Phi_{\delta}^{\epsilon}(s)$ satisfies, for any given $s$ such that $0<s<1$ properties (2.6) and (2.7). Indeed, substituting into (2.7) the value function (2.12) shows that property holds provided

[TABLE]

Inequality (5.26) holds since, for any given $0<s<1$, the function

[TABLE]

is a convex function of $\delta$ that is equal to zero at the boundaries of the interval $0\leq\delta\leq 1$, and this clearly implies $Z_{s}(\delta)\leq 0$ on this interval, independently of the value of $s$.

Property (2.6) holds since the value function $\Phi_{\delta}^{\epsilon}(s)$ is concave with respect to $s\geq 0$. Direct computations allow to conclude that the sign of the second derivative of the value function is non positive if

[TABLE]

Inequality (5.27) holds considering that, for $x\geq 0$, the function

[TABLE]

is a convex function, positive in the interval $0<x\leq(1-\delta)/\epsilon$.

6. Conclusions

We introduced and discussed in the present paper kinetic models able to describe the statistical distribution of alcohol consumption in a multi-agent society. We succeeded to explain the emerging of a generalized Gamma distribution by resorting to the value function theory pioneered by Kahneman and Twersky [41, 42]. In all cases, this macroscopic behavior is a consequence of the choice made at the microscopic level, choice that takes into account the essential features of the human behavior related to the phenomenon of alcohol consumption. The kinetic modeling is similar to the one introduced in [38], subsequently generalized in [39], in which the human behavior was responsible of the formation of a macroscopic equilibrium in the form of a Lognormal distribution. Hence, the present results can be considered as an extension of the kinetic description of [39], which allows to classify, at a microscopic level, the main differences in the elementary interaction which produce a whole class of generalized Gamma distributions, ranging from the classical Gamma density to the Lognormal one. Well-known arguments of kinetic theory allow to model these phenomena by means of a Fokker–Planck equation with variable coefficients of diffusion and drift. Numerical experiments show the perfect agreement between the steady states of the kinetic model of Boltzmann type, and its continuous counterpart, represented by a Fokker–Planck type equation.

We hope that the present analysis about the formation of generalized Gamma distribution in alcohol consumption can be of help to previous studies, like the ones presented in [43, 57, 59], whose analysis was limited to Lognormal, Gamma and Weibull distributions. As remarked in [43], the accurate modeling of alcohol consumption as an upshifted distribution will provide public health decision-makers with accurate data to assess the impact of alcohol consumption within and across countries and will aid in determining public health priorities and where to allocate resources.

Acknowledgement

This work has been written within the activities of GNFM group of INdAM (National Institute of High Mathematics), and partially supported by MIUR project “Optimal mass transportation, geometrical and functional inequalities with applications”. The research was partially supported by the Italian Ministry of Education, University and Research (MIUR): Dipartimenti di Eccellenza Program (2018–2022) - Dept. of Mathematics “F. Casorati”, University of Pavia.