2c:null 2d:["$","title",null,{"children":"Maximum likelihood for iid Poisson data | Mattias Villani | Observable"}] 32:T8da,### The Poisson model and likelihood function Since `NBidders` is a count variable, a natural model to try is an iid Poisson model ${tex.block`X_1, X_2, \ldots, X_n \overset{\mathrm{iid}}{\sim} \text{Pois}(\lambda),`} where ${tex`X_i`} is the number of bidders in auction *i*. The **maximum likelihood** (**ML**) method estimates the unknown parameter ${tex`\lambda`} by maximizing the likelihood function. The **likelihood function** is the probability of the observed data considered as function of the parameter: ${tex.block`L(\lambda) = P(x_1,x_2,\ldots,x_n \vert \lambda)=\prod_{i=1}^n P(x_i \vert \lambda),`} where ${tex`P(x_i \vert \lambda)`} is the [Poisson probability mass function](poisson-distribution). The ML estimator makes a lot of sense: a good value for ${tex`\lambda`} should make the observed data as probable as possible. Conversely, a ${tex`\lambda`} that gives a very low probability to the observed data is probably far from the true/population ${tex`\lambda`} that generated the data. After all, the observed data is what "nature" generated (note that this is all under the assumption that nature's data generating process can be accurately modeled by independent draws from a Poisson distribution). It turns out to be more convenient mathematically to maximize the **log-likelihood** function ${tex.block`\log L(\lambda) = \sum_{i=1}^n \log P(x_i \vert \lambda),`} which gives exactly the same estimate for ${tex`\lambda`} since the logarithm is a monotonically increasing function. The likelihood values are actually not so interpretable. This is because ${tex`L(\lambda)`} is the *joint* probability of all observations in the dataset, which naturally gets smaller and smaller as we obtain more data. In the graph below we therefore plot the *average* probability of a single data point. The natural average for probabilities is the [geometric mean](https://en.wikipedia.org/wiki/Geometric_mean) and we plot the **average probability for a data point** in this geometric sense: ${tex.block`\bar L(\lambda) = \prod_{i=1}^n P(x_i \vert \lambda)^{\frac{1}{n}}`} Note that it does not matter if we maximize ${tex`L(\lambda)`}, ${tex`\log L(\lambda)`} or ${tex`\bar L(\lambda)`}, we will get the same estimate for ${tex`\lambda`}.33:T723,### Interactive exploration of the likelihood function The interactive graph to the left below shows how well a Poisson model with a certain ${tex`\lambda`} fits the observed data (blue bars). The height of the bars are *proportions* of observations falling into that bin, so that they can be compared to the Poisson model's probabilities given by the dots on the orange line. A well fitting Poisson model should therefore have an orange line that is a close fit to the blue bars. With the initial choice of ${tex`\lambda=2`}, we see that the model fits the data poorly. For example, the model gives a probability for zero bidders of around 0.13, but auctions without bidders were observed in less than 1% of the auctions. So the model with ${tex`\lambda=2`} grossly overestimates the probability of zero bidders. Change ${tex`\lambda`} until the model fits the data well. Click on *show ML fitting* to add an optimal model fit according to the maximum likelihood (ML) method. For the Poisson distribution, this ML estimator of ${tex`\lambda`} is actually the mean ${tex`\bar X`} of the data. We prove this mathematically at the end of the notebook. As you can see by interacting with the graph, each value of ${tex`\lambda`} provides an *average probability* for the data; this is the geometric mean ${tex`\bar L(\lambda)`} discussed above. In the figure to the right, we plot all these values as a function of ${tex`\lambda`}. The maximum likelihood estimate of ${tex`\lambda`} is the value that maximizes this function. The orange dot in the graph indicates the value of ${tex`\lambda`} that you have chosen with the slider. When it is at the peak of the function to the right, you have found the maximum likelihood estimate! You can also click the *show log-likelihood* button to instead show the log-likelihood function.34:T5d5,> *This section requires a bit more mathematics. You need to know about the technique from calculus of maximizing a function by finding the point where the derivative of the function is zero.* In the figures above, we maximized the geometric mean of the probability for the observed data. Maximizing this average probability is equivalent to maximizing the joint probability for the entire dataset, ${tex.block`P(x_1,x_2,\ldots,x_n \vert \lambda)=\prod_{i=1}^n P(x_i \vert \lambda)`} where we have used that the observations are assumed to be independent and therefore the joint probability becomes a product of the probability for each observation, ${tex`P(x_i \vert \lambda)`}. Note that data is written with lowercase letters because here we are interested in maximizing the probability for the data material that we have **observed**. In the expression above, we have explicitly written out that the probability of observing a certain data material depends on ${tex`\lambda`}, which is the parameter that we must choose so that the probability for the data is maximized. The joint probability for the data ${tex`P(x_1,x_2,\ldots,x_n \vert \lambda)`} considered as a function of the parameter ${tex`\lambda`} for a fixed (observed) data material is called the **likelihood function**, or just **likelihood**. The maximum likelihood estimate of a model parameter (for example, ${tex`\lambda`} in the Poisson model) is the value of the parameter that maximizes the likelihood function. 35:T50c,In the case where data comes from a Poisson distribution, the probability function for an observation is ${tex`P(x) = \frac{e^{-\lambda}\lambda^x}{x!}`}, and therefore ${tex`\log P(x) = -\lambda + x\log\lambda - \log x!`}, so the log-likelihood is: ${tex.block`\ell(\lambda)=\sum_{i=1}^n\Big(-\lambda + x_i\log\lambda - \log(x_i !) \Big) = -n\lambda + \log\lambda \sum_{i=1}^n x_i - \sum_{i=1}^n \log(x_i !)`} We know from mathematical analysis that the maximum of a function ${tex`f(x)`} is found by setting the first derivative to zero and solving for ${tex`x`}. The first derivative has a simple form: ${tex.block`\frac{d}{d\lambda}\ell(\lambda) = -n + \frac{\sum_{i=1}^n x_i}{\lambda} = 0`} which gives the solution ${tex`\lambda = \frac{\sum_{i=1}^n x_i}{n}=\bar x`}. We can verify that this is indeed a maximum by checking whether the second derivative is negative at ${tex`\lambda = \bar x`}. The second derivative is ${tex.block`\frac{d^2}{d\lambda^2}\ell(\lambda) = - \frac{\sum_{i=1}^n x_i}{\lambda^2},`} which is negative for all ${tex`\lambda`} since both the data and ${tex`\lambda`} must be positive. The maximum likelihood estimator of the parameter ${tex`\lambda`} in the univariate Poisson model is therefore the sample mean ${tex`\hat\lambda = \bar x`}. 😅🥳36:[{"id":"47a98ee79de7b4b3716e827299fc52ebaed0d3bc6bc47841a6eb324fe4c9c19eda2424b7e2b38ecdc5e4c0c332a1112f9df83e2b55f4be3c95dd11f8b5c929e9","href":"https://static.observableusercontent.com/files/47a98ee79de7b4b3716e827299fc52ebaed0d3bc6bc47841a6eb324fe4c9c19eda2424b7e2b38ecdc5e4c0c332a1112f9df83e2b55f4be3c95dd11f8b5c929e9","name":"ebaydata@2.csv","type":"text/csv","lastModified":1707039765018.909,"size":47654}] 30:{"id":"c5e41c4761be764f","slug":"maximum-likelihood-for-iid-poisson-data","title":"Maximum likelihood for iid Poisson data","thumbnail":"4562f7b5839d8333aabf42a2802ef89a976f9745cd0e9d56b0c7cfde6790cba3","access_level":"public","create_time":"$D2024-02-04T09:42:45.018Z","update_time":"$D2025-04-15T13:47:42.573Z","view_time":null,"first_public_time":"$D2024-02-04T11:53:43.763Z","fork":{"id":"630c0c263bea1ca0","slug":"maximum-likelihood-poisson-regression","title":"Maximum likelihood - Poisson regression","owner":{"id":"7592a7ab121f74a0","login":"mattiasvillani","name":"Mattias Villani"}},"stars":0,"starred":false,"owner":{"type":"team","id":"7592a7ab121f74a0","login":"mattiasvillani","name":"Mattias Villani","avatar_url":"https://avatars-next.observableusercontent.com/a440aaf39a3abefb86cbcc758fd7a5ab42555ae68d6825b52ec006f385a94d1d","tier":"pro_2024"},"version":2430,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"56499be8cf43f0ef","login":"mattiasvillani","name":"Mattias Villani","avatar_url":"https://avatars-next.observableusercontent.com/a440aaf39a3abefb86cbcc758fd7a5ab42555ae68d6825b52ec006f385a94d1d","tier":null}],"members":[{"type":"individual","id":"56499be8cf43f0ef","login":"mattiasvillani","name":"Mattias Villani","avatar_url":"https://avatars-next.observableusercontent.com/a440aaf39a3abefb86cbcc758fd7a5ab42555ae68d6825b52ec006f385a94d1d","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"subscribed":null,"forks":1,"collection":{"id":"9fdb7c9ba7be6ba8","slug":"bayesian-learning","title":"Bayesian Learning","owner":{"id":"7592a7ab121f74a0","login":"mattiasvillani","name":"Mattias Villani"}},"type":"notebook","body":{"title":"Maximum likelihood for iid Poisson data","cells":[{"id":38,"value":"# Maximum likelihood for iid Poisson data","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":771,"value":"### Data from coin auctions on eBay\nIn this notebook, we analyze the number of bidders (`NBidders`) in 550 eBay auctions for collector coins^{[[1](#footnote1)]} using the Poisson distribution. The original data contains 1000 auctions, but we will here restrict attention to the auctions with low *reserve prices*, which is the lowest price the seller is willing to accept for the coin. This makes it (more) reasonable to assume that the auctions are independent and identically distributed (iid) from a Poisson distribution with the *same* ${tex`\\lambda`}.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":1158,"value":"$32","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":1175,"value":"$33","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":118,"value":"viewof settings = Inputs.form([\n Inputs.range([0.1, 10], {value: 2, step: 0.1, label: tex`\\lambda:`}),\n Inputs.toggle({label: \"show ML fit\", value: false}),\n Inputs.toggle({label: \"show log-likelihood\", value: false}),\n ])","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":2412,"value":"viewof dataselect = Inputs.radio([\"low reservation prices\",\"high reservation prices\", \"all\"], {label: \"Select auctions:\", value: \"low reservation prices\"})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1121,"value":"textslh = {\n if (!settings[1]){\n return html`${textquantity} for the ${nObs} observed data points with the model \n Pois(λ = ${lambda}) is ${geomean_like(y,lambda).toPrecision(4)}`\n }\n else{\n return html`${textquantity} for the ${nObs} observed data points with the model \n Pois(λ = ${lambda}) is ${geomean_like(y,lambda).toPrecision(4)}
\n ${textquantity} for observed data with the ML estimate \n Pois(λ = ${mle.toPrecision(4)}) is ${geomean_like(y,mle).toPrecision(4)}\n`\n }\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1218,"value":"bothplots = html`\n \n \n \n

\n ${pltdata}\n

\n ${pltlikelihood}\n

\n \n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":2382,"value":"### Poisson regression and beyond\nThe notebook [Poisson regression](maximum-likelihood-poisson-regression) analyzes the dataset with all 1000 auctions using the reserve price as a covariate. Yes, you read correctly, it is possible to analyze count data using regression methods based on the Poisson distribution! It turns out that the maximum likelihood estimator in Poisson regression cannot be derived mathematically, like we do here for the iid model in the next section. But we can obtain the ML estimate in any model quite easily by using numerical optimization software. 🤩","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":1233,"value":"### Derivation of the maximum likelihood estimator","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":1235,"value":"$34","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":1274,"value":"It often turns out to be easier to maximize the *logarithm* of the likelihood function, the so-called **log-likelihood**:\n\n${tex.block`\\ell(\\lambda) = \\log L(\\lambda) = \\log P(x_1,x_2,\\ldots,x_n \\vert \\lambda)=\\sum_{i=1}^n \\log P(x_i \\vert \\lambda)`}\n\nSince the logarithm is a monotonically increasing function, the same \$ \\lambda \$ will maximize both the likelihood and the log-likelihood. The log-likelihood for independent observations is easier to work with because observations enter additively (sum) rather than multiplicatively (product) on the log scale. And it's easier to calculate the derivative of a sum than the derivative of a product.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":1277,"value":"$35","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":1303,"value":"footnote1 = md`1. [^](#paragraph1) Datamaterialet kommer från artikeln: Wegmann, B. och Villani, M. (2011). Bayesian Inference in Structural Second-Price Common Value Auctions, *Journal of Business and Economic Statistics* [pdf](https://www.tandfonline.com/doi/abs/10.1198/jbes.2011.08289)`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":109,"value":"
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n\n---\n#### 🙈 Appendix 🙈\n\nHere is just some code and details that you probably don't want to see.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":221,"value":"#### Imports","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":2047,"value":"import {toc} from \"@mbostock/toc\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":26,"value":"jstat = require('jstat')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":877,"value":"mathjs = require('mathjs')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1300,"value":"footnote = (referencingCellName, number, text) => {\n return html`

^ ${text}

`;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1502,"value":"import {addColumn} from \"@observablehq/excel-add-columns\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1762,"value":"optimjs = require('optimization-js@1.5.0/dist/optimization.js')","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1231,"value":"#### Hidden graphs, shown above","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":780,"value":"pltdata = Plot.plot({\n x: {label: \"number of bidders\"},\n y: {label: \"probability/proportion\"},\n marks: [\n Plot.barY(data, Plot.groupX({y: \"proportion\"}, {x: \"NBidders\", fill: \"steelblue\", opacity: 0.7})),\n Plot.line(poispdf_mle, {filter: d => settings[1], x: \"x\", y: \"pdf\", stroke: \"#5F9EA0\", strokeWidth: 3, \n marker: \"circle\", opacity: 1}),\n Plot.line(poispdf, {x: \"x\", y: \"pdf\", stroke: \"#F28C28\", strokeWidth: 3, marker: \"circle\"}),\n Plot.ruleY([0])\n ]\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1062,"value":"pltlikelihood = Plot.plot({\n y: {label: yaxislabel},\n x: {label: \"λ\"},\n marks: [\n Plot.ruleX([mle], {filter: d => settings[1], stroke: \"#5F9EA0\", strokeWidth: 2}),\n //Plot.ruleY([0]),\n Plot.lineY(meanprob_grid, {x: \"lambda\", y: \"geomean\", stroke: \"steelblue\", strokeWidth: 3}),\n Plot.dot([{x: lambda, y: geomean_like(y, lambda)}], {x: \"x\", y: \"y\", fill: \"orange\", marker: \"circle\", r: 7})\n ]\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":223,"value":"#### Data","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":2374,"value":"reserve_threshold = 0.575 // threshold to obtain the auction with lowest reserve prices","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":2350,"value":"yaxislabel = {if (settings[2]){return \"log-likelihood\"}else{return \"average probability for a single observation\"}}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":2369,"value":"textquantity = {if (settings[2]){return \"Log-likelihood\"}else{return \"Average probability\"}}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1881,"value":"lambda = settings[0]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":768,"value":"ebaydata = FileAttachment(\"ebaydata@2.csv\").csv({typed: true}); ","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":948,"value":"data = {\n if (dataselect == \"low reservation prices\"){\n return ebaydata.filter(d => d.ReservePriceFrac <= reserve_threshold)\n }else{\n if (dataselect == \"high reservation prices\"){\n return ebaydata.filter(d => d.ReservePriceFrac > reserve_threshold)}\n }\n return ebaydata\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":2423,"value":"nObs = data.length","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1900,"value":"data_high = ebaydata.filter(d => d.ReservePriceFrac > reserve_threshold)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":17,"value":"poispdf = {\n const x = d3.range(0, d3.max(ebaydata, d => d.NBidders) + 5, 1);\n const data = x.map(x => ({x: x, pdf: jstat.poisson.pdf(x, lambda)}));\n return data\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":978,"value":"poispdf_mle = {\n const x = d3.range(0, d3.max(y) + 5, 1);\n const data = x.map(x => ({x: x, pdf: jstat.poisson.pdf(x, mle)}));\n return data\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1081,"value":"y = Array.from(d3.map(data, d => d.NBidders))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1904,"value":"y_high = Array.from(d3.map(data_high, d => d.NBidders))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":981,"value":"mle = d3.mean(y)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1902,"value":"mle_high = d3.mean(y_high)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1327,"value":"meanprob_grid = {\n const lambdaGrid = d3.range(1, 10, 0.01);\n const data2 = lambdaGrid.map(lambda => ({lambda: lambda, geomean: geomean_like(y, lambda)}));\n return data2;\n}\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":961,"value":"#### Functions","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":849,"value":"function geomean_like(data, param){\n if (!settings[2]){\n const probs = data.map(x => (Math.exp(-param)*param**x)/mathjs.factorial(x));\n const logprobs = probs.map(d => Math.log(d))\n return Math.exp(mathjs.mean(logprobs))\n }\n else{\n const probs = data.map(x => (Math.exp(-param)*param**x)/mathjs.factorial(x));\n const logprobs = probs.map(d => Math.log(d))\n return d3.sum(logprobs)\n }\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":229,"value":"#### Style","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":202,"value":"","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined"},{"id":58,"value":"viewof md_colour = Inputs.color({label: \"Highlight color cells\", value: \"#F5F5DC\"})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"}],"files":"$W36","stdlib":"1","titling":"auto"}} 31:"$30:body"