32:T931,datasetPure = { function incomeDistribution(nGroups) { let data = [] let groupPopulation = populationTotal/nGroups let groupIncomeMin = 1e-12, groupIncomeMax, groupIncomeMean let groupTotalYearIncome = [] //The function f and its derivative dfdx come in handy for Newtonian iterations let f = function(groupIncomeMin, groupIncomeMax) { return integral(populationWithIncome, groupIncomeMin, groupIncomeMax, 1) - groupPopulation } let dfdx = function(groupIncomeMax) { return populationWithIncome(groupIncomeMax) } let delta0 = 10000 //initial step for (let i = 0; i < nGroups; i++) { if (i < nGroups-1) { // groupIncomeMax = groupIncomeMin + delta0 let groupIncomeMaxPre = groupIncomeMax let err = 100 let j = 0 //Newton iterations: Calculating the maximum income difference at which the population size is equal to a given population size. while (err*err > 0.01) { groupIncomeMax += -f(groupIncomeMin, groupIncomeMaxPre)/dfdx(groupIncomeMaxPre) err = groupIncomeMax - groupIncomeMaxPre groupIncomeMaxPre = groupIncomeMax j++ if (j>15) { groupIncomeMax = j + " iterations were not enough, last result: " + groupIncomeMax groupIncomeMean = j + " iterations were not enough" break } groupTotalYearIncome[i] = integral(incomeOfPopulationGroup, groupIncomeMin, groupIncomeMax, 1)*12 groupIncomeMean = groupTotalYearIncome[i]/12.0/groupPopulation delta0 = groupIncomeMax - groupIncomeMin } } else { //The maximum income can be infinitely large, so we just write that it is very large. groupTotalYearIncome[i] = incomeTotal for (let j = 0; j < nGroups - 1; j++) { groupTotalYearIncome[i] -= groupTotalYearIncome[j] } groupIncomeMean = groupTotalYearIncome[i]/12/groupPopulation groupIncomeMax = 2*groupIncomeMean - groupIncomeMin } data.push({population: groupPopulation, incomeMin: groupIncomeMin, incomeMean: groupIncomeMean, incomeMax: groupIncomeMax, populationIncomeYear: groupTotalYearIncome[i]}) groupIncomeMin = groupIncomeMax } return data } return incomeDistribution }33:[] 30:{"title":"Calculation of income distribution based on the lognormal model","cells":[{"id":1069,"value":"# Calculation of income distribution based on the lognormal model","pinned":false,"hidden":false,"mode":"md"},{"id":1424,"value":"## Inputs","pinned":false,"hidden":false,"mode":"md"},{"id":1184,"value":"Russia, 2020","pinned":false,"hidden":false,"mode":"md"},{"id":1083,"value":"Total population:","pinned":false,"hidden":false,"mode":"md"},{"id":16,"value":"populationTotal=0.5*(146748590 + 146171015)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1087,"value":"Total income (RUB):","pinned":false,"hidden":false,"mode":"md"},{"id":208,"value":"incomeTotal = 62700.8e9","pinned":true,"hidden":false,"mode":"ojs"},{"id":1093,"value":"Mean monthly income per capita, *x_mean* (RUB):","pinned":false,"hidden":false,"mode":"md"},{"id":9,"value":"incomeMean = incomeTotal/populationTotal/12","pinned":true,"hidden":false,"mode":"ojs"},{"id":1102,"value":"Median monthly income per capita, *x_median* (RUB):","pinned":false,"hidden":false,"mode":"md"},{"id":6,"value":"incomeMedian = 27036.4","pinned":true,"hidden":false,"mode":"ojs"},{"id":1428,"value":"## Results","pinned":false,"hidden":false,"mode":"md"},{"id":661,"value":"dataset15 = datasetPure(15)","pinned":true,"hidden":false,"mode":"ojs"},{"id":665,"value":"dataset146 = datasetPure(146)","pinned":true,"hidden":false,"mode":"ojs"},{"id":667,"value":"dataset1465 = datasetPure(1465)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1457,"value":"dataset1465[1464]","pinned":true,"hidden":false,"mode":"ojs"},{"id":1433,"value":"## Theory","pinned":false,"hidden":false,"mode":"md"},{"id":1108,"value":"Lognormal population distribution *f* of some monthly income *x*:","pinned":false,"hidden":false,"mode":"md"},{"id":1111,"value":"f(x) = \\frac1{x\\;\\sigma\\;\\sqrt{2\\;\\pi}}\\;exp\\left(-\\frac{(\\ln x\\;-\\mu\\;)^2}{2\\sigma^2}\\right)","pinned":false,"hidden":false,"mode":"tex"},{"id":0,"value":"populationWithIncome = function(x) {\n return populationTotal/(x*Math.sqrt(2*Math.PI*sigma*sigma))*Math.exp(-Math.pow((Math.log(x)-mu),2)/(2*sigma*sigma))\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1140,"value":"where:","pinned":false,"hidden":false,"mode":"md"},{"id":1127,"value":"\\mu\\;=\\;\\ln(x_{median})","pinned":false,"hidden":false,"mode":"tex"},{"id":12,"value":"mu = Math.log(incomeMedian)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1133,"value":"\\sigma\\;=\\;\\sqrt{2\\;(-\\mu\\;+\\;\\ln(x_{mean}))}","pinned":false,"hidden":false,"mode":"tex"},{"id":14,"value":"sigma = Math.sqrt(2*(-mu + Math.log(incomeMean)))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1216,"value":"fplot = {\n let data = []\n let nPoints = 1000\n let xMin = 1e-12, xMax = 100e3\n for (let i = 0; i < nPoints; i++) {\n let xTemp = xMin + (xMax - xMin)/nPoints*i\n data.push({income: xTemp, population: populationWithIncome(xTemp)})\n }\n return data\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1230,"value":"xMeanLine = {\n let data = []\n data[0] = {income: incomeMean, population: 0}\n data[1] = {income: incomeMean, population: populationWithIncome(incomeMean)}\n return data\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1271,"value":"xMedianLine = {\n let data = []\n data[0] = {income: incomeMedian, population: 0}\n data[1] = {income: incomeMedian, population: populationWithIncome(incomeMedian)}\n return data\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1348,"value":"dotsMeanMedian = {\n let data = []\n data[0] = {income: incomeMean, population: populationWithIncome(incomeMean), label: \"Mean\"}\n data[1] = {income: incomeMedian, population: populationWithIncome(incomeMedian), label: \"Median\"}\n return data\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1214,"value":"Plot.plot({\n x: {\n domain: [0, 101000],\n axis: \"bottom\",\n grid: true,\n label: \"mounthly income (per capita), RUB\",\n labelAnchor: \"right\"\n },\n y: {\n axis: null,\n domain: [0, 4000]\n },\n marks: [\n Plot.line(xMedianLine, {x: \"income\", y: \"population\", stroke: \"#FEAB9A\", text:\"Median\", textAnchor: \"end\"}),\n Plot.line(xMeanLine, {x: \"income\", y: \"population\", stroke: \"#FFD333\"}),\n Plot.line(fplot, {x: \"income\", y: \"population\", stroke: \"#56927E\", strokeWidth: 4}),\n Plot.dot(dotsMeanMedian, {x: \"income\", y: \"population\", r: 5, strokeWidth: 3, stroke: \"#56927E\", fill: \"white\"}),\n Plot.text(dotsMeanMedian, {x: \"income\", y: \"population\", text: \"label\", dy: -10, dx: 12})\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1342,"value":"Let's plot a Lorenz curve based on this data and compare it with Rosstat data.","pinned":false,"hidden":false,"mode":"md"},{"id":1192,"value":"lorenzCurve = {\n let data = []\n let nPoints = 1000\n let xMin = 1e-12, xMax = 1e6\n for (let i = 0; i < nPoints; i++) {\n let xTemp = xMin + (xMax - xMin)/nPoints*i\n data.push({\n population: integral(populationWithIncome, 1e-12, xTemp, 100)/populationTotal,\n income: integral(incomeOfPopulationGroup, 1e-12, xTemp, 1000)*12/incomeTotal\n })\n }\n return data\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1311,"value":"lorenzCurveRosstat = {\n return [\n {population: 0.2, income: 0.055},\n {population: 0.4, income: 0.158},\n {population: 0.6, income: 0.311},\n {population: 0.8, income: 0.538},\n {population: 1, income: 1},\n ]\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1189,"value":"Plot.plot({\n y: {\n grid: true,\n label: \"income\"\n },\n x: {\n grid: true,\n label: \"population\"\n },\n marks: [\n Plot.line(lorenzCurve, {x: \"population\", y: \"income\", stroke: \"#56927E\", strokeWidth: 3}),\n Plot.dot(lorenzCurveRosstat, {x: \"population\", y: \"income\", stroke: \"#56927E\", strokeWidth: 3, fill: \"white\", r: 5})\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":202,"value":"incomeOfPopulationGroup = function(x) {\n return x*populationWithIncome(x)\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1442,"value":"## Dataset generation","pinned":false,"hidden":false,"mode":"md"},{"id":1446,"value":"The task is to determine what is the minimum, mean and maximum income for a group of a given size. The size of the group is given by the number of groups *nGroups* into which the total population is divided.","pinned":false,"hidden":false,"mode":"md"},{"id":326,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"},{"id":90,"value":"## Additional Functions","pinned":false,"hidden":false,"mode":"md"},{"id":1403,"value":"Numerical integration from *a* to *b* with accuracy *tol* by the trapezoidal rule","pinned":false,"hidden":false,"mode":"md"},{"id":45,"value":"integral = function(f, a, b, tol) {\n let n = 1\n let h = b - a\n let sum = []\n let i = 1;\n let s, result\n sum[0] = 0.5*h*(f(a) + f(b))\n let err = 10*tol\n while ((err > tol)) {\n h = h/2\n s = 0\n for (let j = 1; j <= n; j++) {\n s += f(a + (2*j-1)*h)\n }\n n = 2*n\n sum[i] = 0.5*sum[i-1] + h*s\n err = Math.abs(sum[i] - sum[i-1])\n result = sum[i]\n i = i+1\n if (i>20) {\n result = \"integration error\"\n break\n }\n }\n return result\n}","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W33","stdlib":"1","titling":"auto"}