In my Econometrics class yesterday, our teacher discussed a sample dataset that measured the amount of money spent per patient on doctor's visits in a year. This excluded hospital visits and the cost of drugs. The context was a discussion of generalized linear models, and for the purposes of Stata, he found that a gamma distribution worked best. Good enough.

But I started thinking about the dataset itself, and the pdf curve "total dollars per patient" might follow. I'm thinking of it as two variables:

1) The number of visits per person (k). Let's use a Poisson distribution with λ=1, so p(0)=0.368, p(1)=0.368, p(2)=0.184, etc.

2) The cost of each visit (x). Let's use a normal curve with a mean of 50 dollars and a SD of 5 dollars. A person's total expenditure would be the sum of normally distributed variables, or x = N(k50, k25) for k=0 to infinity visits. So a person with no visits would spend nothing, ~.95 of people with one visit would spend between 40 and 60, ~.95 of those with two visits would spend between ~85.86 and ~114.14, three visits would spend ~132.68 and ~167.32, etc.

I think a pdf curve would be some combination of these two functions. My guess is that it would look like a vertical line a 0 and a series of decreasing, ever flattening humps at 50, 100, 150 dollars, etc, where the integral of x=(0, inf.) equals 0.632 (because zero visits is p=.0368)

Can anyone help me put the pieces together? What would that function be? I've taken through Calc III (multivariate/vector calculus), FYI. Thanks.

  • $\begingroup$ I think you are asking the following: Let $\{X_i\}_{i=1}^{\infty}$ be i.i.d. Gaussian $N(m,\sigma^2)$ (with $m=50, \sigma^2=25$). Let $K$ be independent of the $\{X_i\}$ variables and Poisson distributed with parameter $\lambda$ (with $\lambda=1$). Define $Z=\sum_{i=1}^K X_i$. Find the CDF and/or PDF of $Z$. Is this your question? $\endgroup$ – Michael Apr 5 at 23:10
  • $\begingroup$ I'll use MathJax next time. Apologies. $\endgroup$ – Matt Mulligan Apr 6 at 0:39
  • $\begingroup$ I think that's close? I'm not quite sure about the subscripts and superscripts. Let's try this: Let $\{X_i\}_{i=0}^{\infty}$ be i.i.d. Gaussian $N(im,i\sigma^2)$ (with $m=50,σ2=25)$. Let K be independent of the {Xi} variables and Poisson distributed with parameter λ (with λ=1). Define $Z=\sum_{i=0}^K X_i$ I changed it slightly because a) the number of visits starts at zero and b) I need to include something recognizing that Xi is the sum of normally distributed random variables, which is $N(m_1+m_2...m_\infty,\sigma_1^2 + sigma_2^2...sigma_\infty^2)$ I don't know if this works, though. $\endgroup$ – Matt Mulligan Apr 6 at 1:39
  • $\begingroup$ @Michael. Because I forgot to tag you in the first post. I'm using the sum of normally distributed random variables because I'm looking for the total a person would spend in a year with i visits. Thanks for answering. $\endgroup$ – Matt Mulligan Apr 6 at 1:51
  • 1
    $\begingroup$ Anyway, in general, you can use the "law of total probability": For any $z \in \mathbb{R}$ we have $$P[Z \leq z] = \sum_{k=0}^{\infty} P[Z\leq z | K=k]P[K=k]$$ and differentiating with respect to $z$ gives $$f_Z(z) = \sum_{k=0}^{\infty} f_{Z|K=k}(z)P[K=k]$$ where in this case it looks like $f_{Z|K=0}(z)=\delta(z)$ (impulse function). $\endgroup$ – Michael Apr 6 at 2:11

There is a minor problem with using the normal distribution here as it has a positive probability of giving negative values

That being said, your idea of what the density in your model might look like is reasonable. Try the following R function:

plotdensity <- function(meanvisitcost, sdvisitcost, meank){
  maxkconsidered <- max(1, qpois(1 - 10^-7, meank))
  kconsidered <- 0:maxkconsidered
  dkconsidered <- dpois(kconsidered, meank)
  dcostgivenk <- function(x, k){ dnorm(x, meanvisitcost*k, sdvisitcost*sqrt(k)) }
  dcost <- function(x){ outer(x, kconsidered, dcostgivenk) %*% dkconsidered }
  curve(dcost, from=0, to=6*meanvisitcost*max(1,meank), n=100001, col="red")
  segments(0, 0, 0, 10^7, col="red")

Then with your parameters of $50,5,1$ you get the following density, where the spike at the left corresponds to a probability (not density) of $e^{-1}\approx 0.368$

plotdensity(50, 5, 1)

density graph plotdensity(50, 5, 1)

while with BruceET's parameters of $100,20,5$ you get the following, where the spike at the left corresponds to a probability (not density) of $e^{-4}\approx 0.018$

plotdensity(100, 20, 5)

enter image description here


I agree with @Michael that you might be thinking about a random sum of random variables. The number of visits $N \sim \mathsf{Pois}(\lambda)$ and the amount spent on the $i$th visit is $X_i \sim \mathsf{Norm}(\mu, \sigma).$

There are formulas for $E(T)$ and $Var(T),$ where $T = \sum_{i=1}^N X_i.$ For the mean it's $E(T) = E(N)E(X) = \lambda\mu.$ There is a more complicated two-term formula for $Var(T),$ sometimes called the 'total variation'. I'll let you look for that in your text, notes, or online.

Simulating (in R statistical software) for a million hypothetical people, with $\lambda =5$, $\mu = 100,$ $\sigma=20,$ we have the following results. (With $m = 10^6$ iterations, you can expect 2 or 3 "significant" digits of accuracy.)

set.seed(2019)  # for reproducibility
m = 10^6;  lam = 5;  mu = 100;  sg = 20
t = numeric(m) 
for(i in 1:m) {
  n = rpois(1, lam);  t[i] = sum(rnorm(n, mu, sg))

mean(t);  var(t);  sd(t);  mean(t > 700)
[1] 500.0381  # aprx E(T) = 500
[1] 51976.34  # aprx V(T)
[1] 227.9832  # aprx SD(T)
[1] 0.185734  # aprs P(T > 700)

hist(t, prob=T, col="skyblue2", main = "Simulated Dist'n of T")

enter image description here

A normal approximation might not the best way to evaluate probabilities such as $P(T > 700).$

Addendum: Revised histogram as per Comment by @r.e.s.:

enter image description here

Also R's default kernel density estimator added to the earlier histogram from $m = 10^6$ iterations, gives a much clearer picture then just the histogram (although, naturally enough, it misses the spike at $0.)$ Additional code following the hist line is shown below:

lines(density(t), col="red", lwd=2)

enter image description here

  • 1
    $\begingroup$ If we put m = 10^7 and use breaks=400 in the hist command, the resulting histogram is practically indistinguishable from @Henry's plot. $\endgroup$ – r.e.s. Apr 7 at 2:43
  • $\begingroup$ @r.e.s.: Thanks for that news. See Addemdum. $\endgroup$ – BruceET Apr 7 at 5:19

(Expanding on a comment by @Michael, this might help some readers understand the R computations in the accepted answer by @Henry.)

Let $$Y= \begin{cases} 0 &\text{ if }K=0\\ X_1 + \ldots + X_K&\text{ if }K>0 \end{cases}$$ where the $X_i$ are i.i.d. Normal random variables with mean $\mu$ and variance $\sigma^2,$ and $K$ is a Poisson random variable.

Then the CDF of $Y$ is $$\begin{align} P(Y\le y)&=\sum_{k\ge 0}^\infty\, P(Y\le y\mid K=k)\,P(K=k)\\ &=P(Y\le y\mid K=0)\cdot P(K=0)\ +\ \sum_{k\gt 0}^\infty\, P(Y\le y\mid K=k)\, P(K=k)\\ &=[0\le y]\cdot P(K=0)\ +\ \sum_{k\gt 0}^\infty\, P(X_1+\ldots+X_k\le y)\, P(K=k)\\ &=[0\le y]\cdot P(K=0)\ +\ \sum_{k\gt 0}^\infty\, \Phi_{}\left({y-k\mu\over \sigma\sqrt{k}}\right)\, P(K=k)\\ \end{align}$$ where [...] are Iverson brackets and $\Phi$ is the Standard Normal CDF; hence, the density function of $Y$ is

$$\begin{align} f_{Y}(y)&=\delta(y)\cdot P(K=0)\ +\ \sum_{k\gt 0}^\infty\,\underbrace{\phi\left({y-k\mu\over \sigma\sqrt{k}}\right)\,{1\over\sigma\sqrt{k}} }_{\text{density of Normal$(k\mu,k\sigma^2)$}}\,P(K=k)\\ \end{align}$$ where $\delta()$ is the Dirac delta function and $\phi$ is the Standard Normal density function.

That is, the density function of $Y$ is the Poisson-weighted sum of a delta function at $y=0$ together with infinitely many Normal densities with increasing variances $1\sigma^2,2\sigma^2,3\sigma^2,...,$ centered respectively at $y=1\mu,2\mu,3\mu,...$


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.