MLEConvexDensityFunction

import%20marimo%0A%0A__generated_with%20%3D%20%220.10.19%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20from%20mosek.fusion%20import%20Model%2C%20Domain%2C%20Expr%2C%20ObjectiveSense%2C%20Var%0A%20%20%20%20import%20mosek.fusion.pythonic%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20np.random.seed(5)%0A%20%20%20%20import%20sys%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20import%20time%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20Domain%2C%0A%20%20%20%20%20%20%20%20Expr%2C%0A%20%20%20%20%20%20%20%20Model%2C%0A%20%20%20%20%20%20%20%20ObjectiveSense%2C%0A%20%20%20%20%20%20%20%20Var%2C%0A%20%20%20%20%20%20%20%20mo%2C%0A%20%20%20%20%20%20%20%20mosek%2C%0A%20%20%20%20%20%20%20%20np%2C%0A%20%20%20%20%20%20%20%20plt%2C%0A%20%20%20%20%20%20%20%20sys%2C%0A%20%20%20%20%20%20%20%20time%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%20%20%20%20This%20notebook%20is%20based%20on%20the%20work%20of%20Terlaky%20and%20Vial%20in%20%5B%22Computing%20Maximum%20Likelihood%20Estimators%20of%20Convex%20Density%20Functions%22%5D(https%3A%2F%2Fdoi.org%2F10.1137%2FS1064827595286578).%20The%20following%20section%20describes%20the%20problem%20definition%2C%20which%20is%20taken%20nearly%20verbatim%20from%20the%20paper.%20You%20can%20also%20read%20about%20this%20problem%20in%20the%20%5BMOSEK%20Modeling%20Cookbook%5D(https%3A%2F%2Fdocs.mosek.com%2Fmodeling-cookbook%2Fpowo.html%23maximum-likelihood-estimator-of-a-convex-density-function).%20We%20also%20explain%20step%20by%20step%20how%20the%20nonlinear%20model%20is%20transformed%20into%20a%20conic%20model%20which%20can%20be%20fed%20into%20MOSEK.%0A%0A%20%20%20%20%20%20%20%20The%20problem%20addressed%20in%20the%20paper%20is%20to%20estimate%20a%20density%20function%20that%20is%20known%20to%20be%20convex.%20Indeed%2C%20the%20paper%20suggests%20using%20a%20maximum%20likelihood%20estimator%2C%20which%20is%20a%20solution%20to%20a%20linearly%20constrained%20optimization%20problem.%20Let%20%5C(%20Y%20%5C)%20be%20the%20real-valued%20random%20variable%20with%20unknown%2C%20convex%20density%20function%20%5C(%20g%20%5C).%20We%20want%20to%20estimate%20%5C(%20g%20%3A%20%5Cmathbb%7BR%7D_%2B%20%5Cto%20%5Cmathbb%7BR%7D_%2B%20%5C)%20from%20observed%20samples.%20%0A%0A%20%20%20%20%20%20%20%20Let%20%5C(%20%5C%7By_1%2C%20%5Cdots%2C%20y_n%5C%7D%20%5C)%20be%20an%20ordered%20sample%20of%20%5C(%20n%20%5C)%20outcomes%20of%20%5C(%20Y%20%5C).%20We%20shall%20assume%20that%20%5C(%20y_1%20%3C%20y_2%20%3C%20%5Cdots%20%3C%20y_n%20%5C).%20The%20estimator%20of%20%5C(%20g%20%5Cgeq%200%20%5C)%20is%20a%20piecewise%20linear%20function%20%5C(%20%5Ctilde%7Bg%7D%20%3A%20%5By_1%2C%20y_n%5D%20%5Cto%20%5Cmathbb%7BR%7D_%2B%20%5C)%20with%20break%20points%20at%20%5C(%20(y_i%2C%20%5Ctilde%7Bg%7D(y_i))%20%5C)%2C%20%5C(%20i%20%3D%201%2C%20%5Cdots%2C%20n%20%5C).%20%0A%0A%20%20%20%20%20%20%20%20Let%20%5C(%20x_i%20%3E%200%20%5C)%2C%20%5C(%20i%20%3D%201%2C%20%5Cdots%2C%20n%20%5C)%2C%20be%20the%20estimator%20of%20%5C(%20g(y_i)%20%5C).%20The%20objective%20is%20to%20maximize%20the%20likelihood%20function%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%5Ctext%7Bmaximize%7D%20%5Cquad%20f(x)%20%3D%20%5Cprod_%7Bi%3D1%7D%5E%7Bn%7D%20x_i.%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20To%20match%20the%20convexity%20requirement%20for%20the%20one-dimensional%20density%20function%2C%20we%20add%20the%20constraint%20that%20the%20slope%20of%20%5C(%20%5Ctilde%7Bg%7D%20%5C)%20is%20nondecreasing.%20This%20is%20written%20as%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%5Cfrac%7B%5CDelta%20x_i%7D%7B%5CDelta%20y_i%7D%20%5Cleq%20%5Cfrac%7B%5CDelta%20x_%7Bi%2B1%7D%7D%7B%5CDelta%20y_%7Bi%2B1%7D%7D%2C%20%5Cquad%20i%20%3D%202%2C%20%5Cdots%2C%20n-1%2C%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20where%20%5C(%20%5CDelta%20x_i%20%3D%20x_i%20-%20x_%7Bi-1%7D%20%5C)%20and%20%5C(%20%5CDelta%20y_i%20%3D%20y_i%20-%20y_%7Bi-1%7D%20%5C).%20%0A%0A%20%20%20%20%20%20%20%20This%20can%20be%20transformed%20into%20the%20following%20constraint%2C%20which%20will%20be%20the%20first%20constraint%20of%20the%20model.%20%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20-%20%5CDelta%20y_%7Bi%2B1%7D%20x_%7Bi-1%7D%20%2B%20(%5CDelta%20y_i%20%2B%20%5CDelta%20y_%7Bi%2B1%7D)%20x_i%20-%20%5CDelta%20y_i%20x_%7Bi%2B1%7D%20%5Cleq%200%2C%20%5Cquad%20i%20%3D%202%2C%20%5Cdots%2C%20n-1%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%0A%20%20%20%20%20%20%20%20We%20also%20have%20the%20requirement%20that%20%5C(%20%5Ctilde%7Bg%7D%20%5C)%20is%20a%20density%20function%3A%20the%20area%20below%20%5C(%20%5Ctilde%7Bg%7D%20%5C)%20must%20sum%20up%20to%20one%2C%20which%20is%20the%20second%20constraint%20of%20the%20model.%20%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%5Csum_%7Bi%3D1%7D%5E%7Bn-1%7D%20%5CDelta%20y_%7Bi%2B1%7D%20%5Cleft(%20%5Cfrac%7Bx_i%20%2B%20x_%7Bi%2B1%7D%7D%7B2%7D%20%5Cright)%20%3D%201.%0A%20%20%20%20%20%20%20%20%5C%5D%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22We%20will%20use%20the%20following%204%20distributions%20to%20test%20our%20implementation.%20In%20the%20applied%20examples%2C%20the%20exponential%20distribution%20will%20be%20used%20to%20generate%20random%20data.%20Other%20distributions%20can%20also%20be%20tested%20easily%20using%20the%20following%20cells%20if%20desired.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20%23%20Define%20density%20functions%2C%20g%2C%20to%20sample%20y(i)%20points%0A%20%20%20%20%23%20The%20density%20is%20e%5E(-z)%0A%20%20%20%20def%20sample_exponential(n)%3A%0A%20%20%20%20%20%20%20%20return%20np.random.exponential(scale%3D1.0%2C%20size%3Dn)%0A%0A%20%20%20%20%23%20Arcsine%20distribution%3A%20range%20(0%20%3C%20z%20%3C%201)%2C%20density%3A%201%20%2F%20(%CF%80%E2%88%9A(z(1-z)))%0A%20%20%20%20def%20sample_arcsine(n)%3A%0A%20%20%20%20%20%20%20%20U%20%3D%20np.random.uniform(0%2C%201%2C%20n)%20%23%20Generate%20uniform%20samples%0A%20%20%20%20%20%20%20%20return%20np.arcsin(U)%20**%202%20%0A%0A%20%20%20%20%23%20Quadratic%20distribution%3A%20range%20(0%20%E2%89%A4%20z%20%E2%89%A4%202)%2C%20density%3A%201%20-%20z%20%2F%202%0A%20%20%20%20def%20sample_quadratic(n)%3A%0A%20%20%20%20%20%20%20%20U%20%3D%20np.random.uniform(0%2C%201%2C%20n)%20%23%20Generate%20uniform%20samples%0A%20%20%20%20%20%20%20%20return%202%20-%20np.sqrt(4%20-%204%20*%20U)%20%20%0A%0A%20%20%20%20%23%20Inverse%20distribution%3A%20range%20(1%20%E2%89%A4%20z)%2C%20density%3A%201%20%2F%20z%C2%B2%0A%20%20%20%20def%20sample_inverse(n)%3A%0A%20%20%20%20%20%20%20%20u%20%3D%20np.random.uniform(0%2C%201%2C%20n)%20%20%23%20Generate%20uniform%20samples%0A%20%20%20%20%20%20%20%20return%201%20%2F%20(1%20-%20u)%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20sample_arcsine%2C%0A%20%20%20%20%20%20%20%20sample_exponential%2C%0A%20%20%20%20%20%20%20%20sample_inverse%2C%0A%20%20%20%20%20%20%20%20sample_quadratic%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20NodeSize%20%3D%20mo.ui.number(40)%0A%20%20%20%20mo.md(f%22Choose%20the%20amount%20you%20want%20to%20generate%20%3A%20%7BNodeSize%7D)%22)%0A%20%20%20%20return%20(NodeSize%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22The%20%24y%24%20values%20are%20assumed%20to%20be%20non-decreasing%2C%20thus%20the%20generated%20values%20are%20sorted.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20NodeSize%2C%0A%20%20%20%20np%2C%0A%20%20%20%20sample_arcsine%2C%0A%20%20%20%20sample_exponential%2C%0A%20%20%20%20sample_inverse%2C%0A%20%20%20%20sample_quadratic%2C%0A)%3A%0A%20%20%20%20%23%20Number%20of%20samples%0A%20%20%20%20n%20%3D%20NodeSize.value%0A%0A%20%20%20%20%23%20Generate%20samples%20from%20different%20density%20functions%0A%20%20%20%20exponential_samples%20%3D%20sample_exponential(n)%0A%20%20%20%20arcsine_samples%20%3D%20sample_arcsine(n)%0A%20%20%20%20quadratic_samples%20%3D%20sample_quadratic(n)%0A%20%20%20%20inverse_samples%20%3D%20sample_inverse(n)%0A%0A%20%20%20%20%23%20Select%20which%20sample%20to%20use%0A%20%20%20%20y%20%3D%20exponential_samples%0A%0A%20%20%20%20%23%20In%20the%20paper%20the%20y's%20are%20assumed%20to%20be%20in%20non-decreasing%20order%0A%20%20%20%20y.sort()%0A%0A%20%20%20%20%23%20As%20parameter%20of%20the%20model%2C%20we%20need%20the%20delta%20values%20array%20where%20delta_y%20%3D%20y%5Bi%5D%20-%20y%5Bi-1%5D%0A%20%20%20%20delta_y%20%3D%20np.diff(y)%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20arcsine_samples%2C%0A%20%20%20%20%20%20%20%20delta_y%2C%0A%20%20%20%20%20%20%20%20exponential_samples%2C%0A%20%20%20%20%20%20%20%20inverse_samples%2C%0A%20%20%20%20%20%20%20%20n%2C%0A%20%20%20%20%20%20%20%20quadratic_samples%2C%0A%20%20%20%20%20%20%20%20y%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%23%20Exponential%20Cone%20Implementation%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%20%20%20%20To%20solve%20the%20described%20problem%2C%20one%20approach%20could%20be%20to%20apply%20an%20exponential%20cone%20transformation.%20To%20achieve%20this%2C%20we%20can%20take%20the%20natural%20logarithm%20of%20the%20objective%20function.%20Although%20this%20addition%20will%20change%20the%20obtained%20objective%20function%20value%2C%20the%20optimal%20solution%20will%20remain%20the%20same.%20Then%2C%20we%20can%20derive%20the%20following%3A%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%5Cln%20%5Cleft(%20%5Cprod_%7Bi%3D1%7D%5E%7Bn%7D%20x_i%20%5Cright)%20%3D%20%20%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%20%5Cln%20(x_i)%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20Then%20our%20objective%20function%20will%20turn%20into%3A%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%5Ctext%7Bmaximize%7D%20%5Cquad%20f(x)%20%3D%20%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%20%5Cln%20(x_i)%20%5C%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%5Ctext%7Bor%20equivalently%7D%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%5Ctext%7Bminimize%7D%20f(x)%20%3D%20%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%20-%5Cln%20(x_i)%0A%20%20%20%20%20%20%20%20%5C%5D%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%20%20%20%20Then%20as%20a%20second%20step%2C%20we%20can%20turn%20the%20model%20into%20the%20following%3A%20%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%5Ctext%7Bminimize%7D%20%5Cquad%20f(x)%20%3D%20%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%20u_i%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20where%20u%20satisfies%3A%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20u_i%20%5Cgeq%20-%5Cln%20x_i%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20The%20rest%20of%20the%20constraints%20in%20the%20original%20model%20remain%20the%20same.%20Now%2C%20we%20have%20an%20additional%20constraint%20for%20u%20definition%20with%20the%20defined%20objective%20function.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%20%20%20%20Finally%2C%20as%20the%20last%20step%20of%20our%20transformation%2C%20we%20model%20the%20logarithm%20with%20an%20exponential%20cone.%20%0A%0A%20%20%20%20%20%20%20%20The%20exponential%20cone%20is%20the%20set%20of%20triples%20%5C(%20(x_1%2Cx_2%2Cx_3)%20%5C)%20satisfying%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20x_1%20%5Cgeq%20x_2%20e%5E%7Bx_3%20%2F%20x_2%7D%2C%20%5Cquad%20x_2%20%5Cgeq%200%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20With%20this%20definition%20on%20hand%2C%20the%20transformation%20can%20easily%20be%20done%20with%20the%20following%20vector%3A%20%0A%0A%20%20%20%20%20%20%20%20%24%24%5Cbegin%7Bbmatrix%7D%0A%20%20%20%20%20%20%20%20x_1%20%5C%5C%0A%20%20%20%20%20%20%20%20x_2%20%5C%5C%0A%20%20%20%20%20%20%20%20x_3%0A%20%20%20%20%20%20%20%20%5Cend%7Bbmatrix%7D%0A%20%20%20%20%20%20%20%20%3D%0A%20%20%20%20%20%20%20%20%5Cbegin%7Bbmatrix%7D%0A%20%20%20%20%20%20%20%20x%20%5C%5C%0A%20%20%20%20%20%20%20%201%20%5C%5C%0A%20%20%20%20%20%20%20%20-u%0A%20%20%20%20%20%20%20%20%5Cend%7Bbmatrix%7D%0A%20%20%20%20%20%20%20%20%24%24%0A%0A%20%20%20%20%20%20%20%20The%20exponential%20cone%20condition%20in%20this%20case%20evaluates%20to%3A%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20(1)%20%5Cquad%20%20x%20%5Cgeq%20e%5E%7B-u%7D%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20(2)%20%5Cquad%20%20%5Cln%20x%20%5Cgeq%20-u%20%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20(3)%20%5Cquad%20%20-%5Cln%20x%20%5Cleq%20u%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20and%20we%20get%20the%20desired%20relationship%20between%20u%20and%20x.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%20%20%20%20We%20end%20up%20with%20the%20following%20final%20model%3A%20%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%5Ctext%7Bmin%7D%20%5Cquad%20f(x)%20%3D%20%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%20u_i%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20-%20%5CDelta%20y_%7Bi%2B1%7D%20x_%7Bi-1%7D%20%2B%20(%5CDelta%20y_i%20%2B%20%5CDelta%20y_%7Bi%2B1%7D)%20x_i%20-%20%5CDelta%20y_i%20x_%7Bi%2B1%7D%20%5Cleq%200%2C%20%5Cquad%20i%20%3D%202%2C%20%5Cdots%2C%20n-1%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%5Csum_%7Bi%3D1%7D%5E%7Bn-1%7D%20%5CDelta%20y_%7Bi%2B1%7D%20%5Cleft(%20%5Cfrac%7Bx_i%20%2B%20x_%7Bi%2B1%7D%7D%7B2%7D%20%5Cright)%20%3D%201.%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%0A%20%20%20%20%20%20%20%20%24%24%0A%20%20%20%20%20%20%20%20%5Cbegin%7Bbmatrix%7D%0A%20%20%20%20%20%20%20%20x%20%5C%5C%0A%20%20%20%20%20%20%20%201%20%5C%5C%0A%20%20%20%20%20%20%20%20-u%0A%20%20%20%20%20%20%20%20%5Cend%7Bbmatrix%7D%20%5Cin%20K_%7BEXP%7D%0A%20%20%20%20%20%20%20%20%24%24%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(Domain%2C%20Expr%2C%20Model%2C%20ObjectiveSense%2C%20time)%3A%0A%20%20%20%20def%20ExponentialCone(delta_y%2Cn%2Cweight)%3A%0A%20%20%20%20%20%20%20%20start%20%3D%20time.time()%0A%20%20%20%20%20%20%20%20M%20%3D%20Model()%0A%0A%20%20%20%20%20%20%20%20%23%20Decision%20Variable%20for%20Conic%20Transformation%20%0A%20%20%20%20%20%20%20%20u%20%3D%20M.variable(%22u%22%2C%20n%20%2CDomain.unbounded())%0A%20%20%20%20%20%20%20%20%23%20The%20decision%20variable%20x%5Bi%5D%2C%20is%20the%20estimator%20of%20g(y%5Bi%5D)%0A%20%20%20%20%20%20%20%20x%20%3D%20M.variable(%22x%22%2C%20n%20%2CDomain.unbounded())%0A%0A%20%20%20%20%20%20%20%20%23%20The%20slope%20of%20g%20is%20ensured%20to%20be%20non-decreasing%2C%20first%20constraint%20%0A%20%20%20%20%20%20%20%20M.constraint(%22Slope%20Integrity%22%2C%20-%20Expr.mulElm(delta_y%5B1%3An%5D%20%2C%20x%5B0%3An-2%5D)%20%2B%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Expr.mulElm((delta_y%5B0%3An-2%5D%20%2B%20delta_y%5B1%3An%5D)%20%2C%20x%5B1%3An-1%5D)%20-%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Expr.mulElm(delta_y%5B0%3An-2%5D%20%2C%20x%5B2%3An%5D)%20%3C%3D%200%20)%0A%0A%20%20%20%20%20%20%20%20%23%20The%20area%20below%20our%20density%20function%2C%20g%2C%20must%20sum%20up%20to%201%2C%20second%20constraint%0A%20%20%20%20%20%20%20%20M.constraint(%22Probability%20Sum%22%2C%20%20Expr.dot(((x%5B0%3An-1%5D%20%2B%20x%5B1%3An%5D)%20%2F%202%20)%20%2C%20delta_y)%20%3D%3D%201.0)%0A%0A%20%20%20%20%20%20%20%20%23%20Exponential%20Cone%20Transformation%20Definition%0A%20%20%20%20%20%20%20%20M.constraint(%22Cone%20Definition%22%2C%20Expr.hstack(x%2CExpr.ones(n)%2C-u)%2C%20Domain.inPExpCone())%0A%0A%20%20%20%20%20%20%20%20M.objective(%22Objective%20Function%22%2C%20ObjectiveSense.Minimize%2C%20Expr.dot(weight%2Cu))%0A%20%20%20%20%20%20%20%20M.solve()%0A%20%20%20%20%20%20%20%20sol%20%3D%20x.level()%0A%20%20%20%20%20%20%20%20ModelStatus%20%3D%20M.getPrimalSolutionStatus()%0A%20%20%20%20%20%20%20%20print(ModelStatus)%0A%20%20%20%20%20%20%20%20end%20%3D%20time.time()%0A%0A%20%20%20%20%20%20%20%20iterations%20%3D%20M.getSolverIntInfo(%22intpntIter%22)%0A%20%20%20%20%20%20%20%20print(f%22Number%20of%20interior-point%20iterations%3A%20%7Biterations%7D%22)%0A%0A%20%20%20%20%20%20%20%20M.writeTask(%22MLE_ExponentialCone.ptf%22)%0A%20%20%20%20%20%20%20%20print(%22Elapsed%20time%3A%20%22%2C%20round(end%20-%20start%2C4)%2C%22%20s.%22)%0A%20%20%20%20%20%20%20%20return%20sol%0A%20%20%20%20return%20(ExponentialCone%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22Let's%20also%20define%20a%20function%20to%20plot%20the%20fit%20of%20our%20estimator.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(plt)%3A%0A%20%20%20%20def%20PlotFigure(y%2Csol)%3A%0A%20%20%20%20%20%20%20%20plt.figure(figsize%3D(10%2C%206))%0A%20%20%20%20%20%20%20%20plt.plot(y%2C%20sol%2C%20label%3D'Y%20Array'%2C%20color%3D'blue'%2C%20marker%3D'o'%2C%20markersize%3D4%2C%20linewidth%3D1)%0A%20%20%20%20%20%20%20%20plt.title('Plot%20of%20Y%20Array')%0A%20%20%20%20%20%20%20%20plt.xlabel('Index')%0A%20%20%20%20%20%20%20%20plt.ylabel('Y%20Values')%0A%20%20%20%20%20%20%20%20plt.grid(True)%0A%20%20%20%20%20%20%20%20plt.legend()%0A%20%20%20%20%20%20%20%20plt.show()%0A%20%20%20%20return%20(PlotFigure%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22The%20retreived%20result%20of%20the%20model%20is%20reported.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(ExponentialCone%2C%20delta_y%2C%20n%2C%20np)%3A%0A%20%20%20%20sol%20%3D%20ExponentialCone(delta_y%2Cn%2Cnp.ones(n))%0A%20%20%20%20return%20(sol%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%22%22%22As%20seen%20in%20the%20figure%2C%20the%20estimator%20outputs%20a%20curve%20resembling%20the%20exponential%20distribution.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(PlotFigure%2C%20sol%2C%20y)%3A%0A%20%20%20%20PlotFigure(y%2Csol)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%23%20Power%20Cone%20Implementation%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%20%20%20%20Alternatively%2C%20one%20could%20also%20solve%20the%20original%20model%20using%20power%20cone%20reformulation.%20Let's%20consider%20again%20maximizing%20the%20likelihood%20function%3A%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%5Ctext%7Bmaximize%7D%20%5Cquad%20f(x)%20%3D%20%5Cprod_%7Bi%3D1%7D%5E%7Bn%7D%20x_i.%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20We%20can%20use%20geometric%20mean%20transformation%20to%20reach%20the%20power%20cone%20representation.%20The%20definition%20of%20the%20power%20cone%20%5C(%20%5Cmathcal%7BP%7D_%7Bn%7D%5E%7B%5Calpha_1%2C%20%5Cdots%2C%20%5Calpha_m%7D%20%5C)%2C%20where%20%5C(%20%5Calpha_i%20%3E%200%20%5C)%20and%20%5C(%20%5Csum%20%5Calpha_i%20%3D%201%20%5C)%7D%20%20is%20the%20following%3A%0A%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20x_1%5E%7B%5Calpha_1%7D%20%5Ccdots%20x_m%5E%7B%5Calpha_m%7D%20%5Cgeq%20%5Csqrt%7Bx_%7Bm%2B1%7D%5E2%20%2B%20%5Cdots%20%2B%20x_n%5E2%7D%2C%20%5Cquad%20x_1%2C%20%5Cdots%2C%20x_m%20%5Cgeq%200%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20Maximizing%20%5C(%20f(x)%20%5C)%20is%20the%20same%20as%20solving%20the%20problem%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%5Ctext%7Bmaximize%7D%20%5Cquad%20t%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%20%5Cquad%20%5Cprod_%7Bi%3D1%7D%5E%7Bn%7D%20x_i%5E%7B1%20%5Cover%20n%7D%20%5Cgeq%20t%0A%20%20%20%20%20%20%20%20%5C%5D%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%20%20%20%20This%20last%20version%20fits%20the%20power%20cone%20definition%2C%20thus%20the%20transformation%20can%20be%20applied%20directly.%20The%20vector%20defining%20the%20power%20cone%20is%20added%20to%20the%20model.%20At%20last%20the%20full%20model%20will%20look%20as%20the%20following%3A%20%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%5Ctext%7Bmaximize%7D%20%5Cquad%20t%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20-%20%5CDelta%20y_%7Bi%2B1%7D%20x_%7Bi-1%7D%20%2B%20(%5CDelta%20y_i%20%2B%20%5CDelta%20y_%7Bi%2B1%7D)%20x_i%20-%20%5CDelta%20y_i%20x_%7Bi%2B1%7D%20%5Cleq%200%2C%20%5Cquad%20i%20%3D%202%2C%20%5Cdots%2C%20n-1%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20%5C%5B%0A%20%20%20%20%20%20%20%20%5Csum_%7Bi%3D1%7D%5E%7Bn-1%7D%20%5CDelta%20y_%7Bi%2B1%7D%20%5Cleft(%20%5Cfrac%7Bx_i%20%2B%20x_%7Bi%2B1%7D%7D%7B2%7D%20%5Cright)%20%3D%201.%0A%20%20%20%20%20%20%20%20%5C%5D%0A%0A%20%20%20%20%20%20%20%20%24%24%0A%20%20%20%20%20%20%20%20%5Cbegin%7Bbmatrix%7D%0A%20%20%20%20%20%20%20%20x_1%20%5C%5C%0A%20%20%20%20%20%20%20%20x_2%20%5C%5C%0A%20%20%20%20%20%20%20%20%5Cdots%20%5C%5C%0A%20%20%20%20%20%20%20%20x_n%20%5C%5C%20%0A%20%20%20%20%20%20%20%20t%0A%20%20%20%20%20%20%20%20%5Cend%7Bbmatrix%7D%0A%20%20%20%20%20%20%20%20%20%5Cin%20%5Cmathcal%7BP%7D_%7Bn%2B1%7D%5E%7B1%2Fn%2C%5Cldots%2C1%2Fn%7D%20%24%24%0A%0A%20%20%20%20%20%20%20%20This%20special%20instance%20of%20the%20power%20cone%20appears%20also%20in%20MOSEK%20as%20the%20%5Bgeometric%20mean%20cone%5D(https%3A%2F%2Fdocs.mosek.com%2Fmodeling-cookbook%2Fpowo.html%23the-power-cone-s).%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(Domain%2C%20Expr%2C%20Model%2C%20ObjectiveSense%2C%20Var%2C%20n%2C%20time)%3A%0A%20%20%20%20def%20PowerCone(delta_y)%3A%0A%20%20%20%20%20%20%20%20start%20%3D%20time.time()%0A%20%20%20%20%20%20%20%20M%20%3D%20Model()%0A%0A%20%20%20%20%20%20%20%20%23%20Decision%20Variable%20for%20Conic%20Transformation%20%0A%20%20%20%20%20%20%20%20u%20%3D%20M.variable(%22u%22%2C%201%20%2CDomain.unbounded())%0A%20%20%20%20%20%20%20%20%23%20The%20decision%20variable%20x%5Bi%5D%2C%20is%20the%20estimator%20of%20g(y%5Bi%5D)%0A%20%20%20%20%20%20%20%20x%20%3D%20M.variable(%22x%22%2C%20n%20%2CDomain.greaterThan(0.0))%0A%0A%20%20%20%20%20%20%20%20%23%20The%20slope%20of%20g%20is%20ensured%20to%20be%20non-decreasing%20%0A%20%20%20%20%20%20%20%20M.constraint(%22Slope%20Integrity%22%2C%20-%20Expr.mulElm(delta_y%5B1%3An%5D%20%2C%20x%5B0%3An-2%5D)%20%2B%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Expr.mulElm((delta_y%5B0%3An-2%5D%20%2B%20delta_y%5B1%3An%5D)%20%2C%20x%5B1%3An-1%5D)%20-%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Expr.mulElm(delta_y%5B0%3An-2%5D%20%2C%20x%5B2%3An%5D)%20%3C%3D%200)%0A%0A%20%20%20%20%20%20%20%20%23%20The%20area%20below%20our%20density%20function%2C%20g%2C%20must%20sum%20up%20to%201%0A%20%20%20%20%20%20%20%20M.constraint(%22Probability%20Sum%22%2C%20%20Expr.dot(((x%5B0%3An-1%5D%20%2B%20x%5B1%3An%5D)%20%2F%202%20)%20%2C%20delta_y)%20%20%3D%3D%201.0)%0A%0A%20%20%20%20%20%20%20%20%23%20Power%20Cone%20Transformation%20Definition%0A%20%20%20%20%20%20%20%20M.constraint(Var.vstack(x%2C%20u)%2C%20Domain.inPPowerCone(%5B1%2Fn%5D%20*%20n))%0A%0A%20%20%20%20%20%20%20%20M.objective(%22Objective%20Function%22%2C%20ObjectiveSense.Maximize%2C%20u)%0A%0A%20%20%20%20%20%20%20%20M.solve()%0A%0A%20%20%20%20%20%20%20%20sol%20%3D%20x.level()%0A%20%20%20%20%20%20%20%20ModelStatus%20%3D%20M.getPrimalSolutionStatus()%0A%20%20%20%20%20%20%20%20print(ModelStatus)%0A%20%20%20%20%20%20%20%20end%20%3D%20time.time()%0A%0A%20%20%20%20%20%20%20%20M.writeTask(%22MLE_PowerCone.ptf%22)%0A%20%20%20%20%20%20%20%20iterations%20%3D%20M.getSolverIntInfo(%22intpntIter%22)%0A%20%20%20%20%20%20%20%20print(f%22Number%20of%20interior-point%20iterations%3A%20%7Biterations%7D%22)%0A%0A%20%20%20%20%20%20%20%20print(%22Elapsed%20time%3A%20%22%2C%20round(end%20-%20start%2C4)%20%2C%22%20s.%22)%0A%20%20%20%20%20%20%20%20return%20sol%0A%20%20%20%20return%20(PowerCone%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22The%20retreived%20result%20of%20the%20model%20is%20reported.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(PowerCone%2C%20delta_y)%3A%0A%20%20%20%20sol_PC%20%3D%20PowerCone(delta_y)%0A%20%20%20%20return%20(sol_PC%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22Once%20more%2C%20as%20seen%20in%20the%20figure%2C%20the%20estimator%20outputs%20a%20curve%20resembling%20the%20exponential%20distribution.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(PlotFigure%2C%20sol_PC%2C%20y)%3A%0A%20%20%20%20PlotFigure(y%2Csol_PC)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22The%20performance%20of%20two%20methods%20is%20evaluated%20using%20different%20sample%20sizes.%20Each%20sample%20is%20solved%20100%20times%20with%20both%20methods%2C%20and%20the%20total%20elapsed%20time%20and%20interior-point%20iterations%20are%20recorded%20in%20the%20table.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20%3C!DOCTYPE%20html%3E%0A%20%20%20%20%20%20%20%20%3Chtml%20lang%3D%22en%22%3E%0A%20%20%20%20%20%20%20%20%3Chead%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%3Cmeta%20charset%3D%22UTF-8%22%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%3Cmeta%20name%3D%22viewport%22%20content%3D%22width%3Ddevice-width%2C%20initial-scale%3D1.0%22%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%3Ctitle%3EFormatted%20Table%3C%2Ftitle%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%3Cstyle%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20table%20%7B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20border-collapse%3A%20collapse%3B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20width%3A%20100%25%3B%20%2F*%20Make%20it%20fit%20the%20container%20*%2F%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20max-width%3A%20900px%3B%20%2F*%20Set%20a%20reasonable%20max%20width%20*%2F%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20margin%3A%20auto%3B%20%2F*%20Center%20the%20table%20*%2F%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%7D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20th%2C%20td%20%7B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20border%3A%201px%20solid%20black%3B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20padding%3A%208px%3B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20text-align%3A%20center%3B%20%2F*%20Ensure%20text%20is%20properly%20aligned%20*%2F%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%7D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20th%20%7B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20background-color%3A%20%23f2f2f2%3B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%7D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20.table-container%20%7B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20display%3A%20flex%3B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20justify-content%3A%20center%3B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20width%3A%20100%25%3B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%7D%0A%20%20%20%20%20%20%20%20%20%20%20%20%3C%2Fstyle%3E%0A%20%20%20%20%20%20%20%20%3C%2Fhead%3E%0A%20%20%20%20%20%20%20%20%3Cbody%3E%0A%0A%20%20%20%20%20%20%20%20%3Cdiv%20class%3D%22table-container%22%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%3Ctable%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctr%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Cth%20rowspan%3D%222%22%20style%3D%22text-align%3A%20center%3B%22%3ENumber%20of%20Samples%3C%2Fth%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Cth%20colspan%3D%222%22%20style%3D%22text-align%3A%20center%3B%22%3EExponential%20Cone%3C%2Fth%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Cth%20colspan%3D%222%22%20style%3D%22text-align%3A%20center%3B%22%3EPower%20Cone%3C%2Fth%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3C%2Ftr%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctr%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Cth%3EElapsed%20Time%3C%2Fth%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Cth%3EInterior-Point%20Iterations%3C%2Fth%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Cth%3EElapsed%20Time%3C%2Fth%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Cth%3EInterior-Point%20Iterations%3C%2Fth%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3C%2Ftr%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctr%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E10%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E0.51%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E900%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E0.55%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E1000%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3C%2Ftr%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctr%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E100%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E1.60%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E2600%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E2.26%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E1800%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3C%2Ftr%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctr%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E500%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E8.75%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E4480%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E11.71%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E2300%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3C%2Ftr%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctr%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E1000%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E16.89%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E5600%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E30.77%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E3100%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3C%2Ftr%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctr%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E2000%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E109.04%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E10600%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E73.32%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E4000%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3C%2Ftr%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctr%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E10000%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E408.44%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E11600%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E241.27%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3Ctd%3E4200%3C%2Ftd%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3C%2Ftr%3E%0A%20%20%20%20%20%20%20%20%20%20%20%20%3C%2Ftable%3E%0A%20%20%20%20%20%20%20%20%3C%2Fdiv%3E%0A%0A%20%20%20%20%20%20%20%20%3C%2Fbody%3E%0A%20%20%20%20%20%20%20%20%3C%2Fhtml%3E%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%20%20%20%20The%20power%20cone%20demonstrates%20significantly%20better%20performance%20in%20terms%20of%20total%20interior-point%20iterations.%20While%20the%20number%20of%20iterations%20increases%20for%20both%20methods%20as%20the%20sample%20size%20grows%2C%20the%20power%20cone%20shows%20a%20slower%20rate%20of%20increase%20compared%20to%20the%20exponential%20cone.%0A%0A%20%20%20%20%20%20%20%20Due%20to%20the%20way%20the%20power%20cone%20is%20currently%20implemented%20in%20Mosek%2C%20the%20time%20per%20interior-point%20iteration%20in%20the%20power%20case%20is%20higher%20than%20in%20the%20exponential%20cone%20case.%20It%20is%20expected%20that%20in%20a%20future%20version%20of%20Mosek%2C%20the%20time%20per%20iteration%20will%20be%20almost%20identical%20for%20both%20formulations.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%23%20Clustering%20Scheme%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%20%20%20%20As%20the%20scale%20of%20the%20problem%20increases%2C%20it%20becomes%20harder%20to%20solve.%20To%20handle%20this%20situation%2C%20a%20clustering%20scheme%20is%20proposed.%20To%20determine%20a%20good%20enough%20estimator%20for%20a%20distribution%2C%20points%20that%20are%20too%20close%20to%20each%20other%20can%20be%20handled%20as%20one.%20While%20the%20points%20are%20merged%20by%20taking%20their%20averages%2C%20we%20add%20their%20weight%20to%20the%20objective%20function%20to%20reflect%20the%20density%20of%20these%20merged%20points.%20As%20a%20result%2C%20we%20end%20up%20using%20less%20points%20in%20our%20model.%20Therefore%2C%20the%20resulting%20clustered%20y%20points%20get%20resorted%2C%20and%20the%20value%20of%20n%20is%20updated.%20%0A%0A%20%20%20%20%20%20%20%20The%20clusters%20are%20built%20in%20a%20manner%20to%20ensure%20that%20the%20difference%20between%20the%20adjacent%20points%20does%20not%20exceed%20the%20resolution.%20The%20resolution%20can%20be%20tuned%20according%20to%20the%20sample%20size%20and%20the%20distribution%20function.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20def%20ClusterPoints(y)%3A%0A%20%20%20%20%20%20%20%20r%20%3D%202e-4%20%23resolution%0A%20%20%20%20%20%20%20%20clustered_y%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20new_cluster%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20w%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20for%20z%20in%20y%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20not%20new_cluster%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20new_cluster.append(z)%0A%20%20%20%20%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20z%20-%20new_cluster%5B-1%5D%20%3C%3D%20r%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20new_cluster.append(z)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20w.append(len(new_cluster))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20clustered_y.append(np.mean(new_cluster))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20new_cluster%20%3D%20%5Bz%5D%0A%0A%20%20%20%20%20%20%20%20%23%20As%20parameter%20of%20the%20model%2C%20we%20need%20the%20delta%20values%20array%20where%20delta_y%20%3D%20y%5Bi%5D%20-%20y%5Bi-1%5D%0A%20%20%20%20%20%20%20%20clustered_y.sort()%0A%20%20%20%20%20%20%20%20delta_y%20%3D%20np.diff(clustered_y)%0A%20%20%20%20%20%20%20%20n%20%3D%20len(clustered_y)%0A%20%20%20%20%20%20%20%20print(%22New%20sample%20set%20length%20went%20down%20from%20%22%2Clen(y)%2C%20%22%20to%20%22%2C%20n%2C%20%22.%22)%0A%20%20%20%20%20%20%20%20return%20clustered_y%2C%20delta_y%2C%20w%2C%20n%0A%20%20%20%20return%20(ClusterPoints%2C)%0A%0A%0A%40app.cell%0Adef%20_(ClusterPoints%2C%20y)%3A%0A%20%20%20%20clustered_y%2C%20clustered_delta_y%2C%20w%2C%20clustered_n%20%3D%20ClusterPoints(y)%0A%20%20%20%20return%20clustered_delta_y%2C%20clustered_n%2C%20clustered_y%2C%20w%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22We%20solve%20the%20model%20once%20more%20with%20the%20new%20points.%20Please%20note%20that%20the%20weight%20array%20is%20not%20a%20dummy%20array%2C%20an%20array%20of%20ones%20as%20it%20was%20in%20previous%20implemention%2C%20and%20is%20actually%20applied%20to%20the%20objective%20function.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(ExponentialCone%2C%20clustered_delta_y%2C%20clustered_n%2C%20w)%3A%0A%20%20%20%20sol_clustered%20%3D%20ExponentialCone(clustered_delta_y%2C%20clustered_n%2C%20w)%0A%20%20%20%20return%20(sol_clustered%2C)%0A%0A%0A%40app.cell%0Adef%20_(PlotFigure%2C%20clustered_y%2C%20sol_clustered)%3A%0A%20%20%20%20PlotFigure(clustered_y%2C%20sol_clustered)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22For%20a%20more%20meaningful%20implementation%2C%20please%20increase%20the%20number%20of%20points%20to%20a%20significantly%20big%20number.%22%22%22)%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A