tag:blogger.com,1999:blog-73791109607960141702017-10-16T20:43:43.207-04:00Clueless FundatmaA random walk through a subset of things I care about. Science, math, computing, higher education, open source software, economics, food etc.Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.comBlogger708125tag:blogger.com,1999:blog-7379110960796014170.post-37398639984647060462017-10-12T08:37:00.000-04:002017-10-12T08:40:05.360-04:00Introduce Concepts in Historical Order?<div dir="ltr" style="text-align: left;" trbidi="on">Let me confess: I have read very few scientific classics in the original.<br /><br />I haven't read the <a href="https://en.wikipedia.org/wiki/Philosophi%C3%A6_Naturalis_Principia_Mathematica">Principia</a>, the <a href="https://en.wikipedia.org/wiki/On_the_Origin_of_Species">Origin of Species</a>, or the <a href="https://en.wikipedia.org/wiki/Euclid%27s_Elements">Elements</a>.<br /><br />I had not even read Einstein's 1905 classic on <a href="https://en.wikipedia.org/wiki/Brownian_motion">Brownian motion</a>, until a few years ago, even though half of my research is directly or indirectly animated by it.<br /><br />Ever since I saw this amazing <a href="https://sachinashanbhag.blogspot.com/2017/09/complex-numbers-part-deux.html">series on complex numbers</a>, I have been wondering whether presenting the historical progression of ideas might be "better" than the standard textbook introduction. Here are some of my observations.<br /><br />The <b>historical approach</b> (HA) is inherently interesting, because it is about ideas and the <b>people</b> behind them. Stories of humans exploring and pushing boundaries, regardless of domain, are fascinating. These stories often have imperfect people grappling with new ideas, getting confused by their implications, arguing back and forth, improving, and gradually perfecting them over centuries. This happened with classical mechanics, evolution, complex numbers, quantum mechanics, etc.<br /><br />The <b>standard approach</b> (SA), on the other hand, steers away from messy pasts, leaps of intuition that came seemingly from nowhere, the entertaining bickering, and the trials and errors. It trims away the excess fat of distractions, consolidates different viewpoints, and presents a sanitized account of an idea. It is, without question, the quickest and cleanest way to learn a new concept. This is an extremely desirable feature in university courses, which have a mandate to "cover" a set of concepts, often in limited time.<br /><br />Perhaps, a good practical compromise is to start with an example rooted in the historical approach to motivate the topic, and transition to the standard textbook approach to teach the meat of the topic. It might be interesting to conclude once again with a historical perspective, perhaps mixed with a discussion of the current state of art and open questions.</div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/jLT6Hx6b_zI" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/10/introduce-concepts-in-historical-order.htmltag:blogger.com,1999:blog-7379110960796014170.post-8182794377619660842017-10-04T14:34:00.000-04:002017-10-04T14:34:07.364-04:00LaTeX: Figure Captions<div dir="ltr" style="text-align: left;" trbidi="on">A minimal working example.<br /><div class="separator" style="clear: both; text-align: center;"><a href="https://2.bp.blogspot.com/-WWLXRbasKs8/WbLiSvJSjFI/AAAAAAAADoQ/lPt8x-WIJ0McAhTWLb3Hrh6dLmnICqYcgCLcBGAs/s1600/Screenshot%2Bfrom%2B2017-09-08%2B14-29-11.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="717" data-original-width="374" src="https://2.bp.blogspot.com/-WWLXRbasKs8/WbLiSvJSjFI/AAAAAAAADoQ/lPt8x-WIJ0McAhTWLb3Hrh6dLmnICqYcgCLcBGAs/s1600/Screenshot%2Bfrom%2B2017-09-08%2B14-29-11.png" /></a></div><br /><br /><script src="https://gist.github.com/shane5ul/3e0a7bd37c1871939ac837a949793d5a.js"></script> <br /><br /></div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/86XZGOs0BRU" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/10/latex-figure-captions.htmltag:blogger.com,1999:blog-7379110960796014170.post-29158010016493949802017-09-25T14:45:00.000-04:002017-09-25T14:45:10.376-04:00Prony Method<div dir="ltr" style="text-align: left;" trbidi="on">Given N equispaced data-points \(F_i = F(t = i \Delta t)\), where \(i = 0, 1, ..., N-1\), <a href="http://sachinashanbhag.blogspot.com/2017/08/exam-question-on-fitting-sums-of.html">the Prony method</a> can be used to fit a sum of m decaying exponenitals: \[F(t) = \sum_{i=1}^{m} a_i e^{b_i t}. \] The 2m unknowns are \(a_i\) and \(b_i\).<br /><div><br /></div><div>In the Prony method, the number of modes in the exponential (m) is pre-specified. There are other methods, which are more general.</div><br />Here is a python subprogram which implements the Prony method.<br /><script src="https://gist.github.com/shane5ul/8b360bb605baa9e29e5e6ede364f4d7d.js"></script> If you have arrays t and F, it can be called as:<br /><br /><span style="font-family: "courier new" , "courier" , monospace;">a_est, b_est = prony(t, F, m)</span></div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/phcwnJW_BZ0" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/09/prony-method.htmltag:blogger.com,1999:blog-7379110960796014170.post-18444359853313708372017-09-22T11:17:00.000-04:002017-09-22T11:17:07.568-04:00MCMC Samplers Visualization<div dir="ltr" style="text-align: left;" trbidi="on">A really nice <a href="https://chi-feng.github.io/mcmc-demo/app.html#HamiltonianMC,banana">interactive gallery of MCMC</a> samplers<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://2.bp.blogspot.com/--dfx5PsCiGw/Wb6RySOpg5I/AAAAAAAADpY/l_zZmiaD-vYjMpJVL2c4fbQbchqybYxhgCLcBGAs/s1600/Screenshot%2Bfrom%2B2017-09-17%2B11-16-13.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="802" data-original-width="1600" height="200" src="https://2.bp.blogspot.com/--dfx5PsCiGw/Wb6RySOpg5I/AAAAAAAADpY/l_zZmiaD-vYjMpJVL2c4fbQbchqybYxhgCLcBGAs/s400/Screenshot%2Bfrom%2B2017-09-17%2B11-16-13.png" width="400" /></a></div><br />You can choose different algorithms, and target distributions, change method parameters and observe the chain evolve.<br /><br />This might come in handy next semester, when I teach a Monte Carlo class.</div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/iC1ZS30g9pA" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/09/mcmc-samplers-visualization.htmltag:blogger.com,1999:blog-7379110960796014170.post-14982650959493721682017-09-19T11:14:00.000-04:002017-09-19T11:14:01.919-04:00Some Useful Math Links!<div dir="ltr" style="text-align: left;" trbidi="on">1. The history of the division symbol (obelus) is fascinating! <a href="https://divisbyzero.com/2017/09/15/the-division-symbol-goes-viral/">(DivisionByZero</a>)<br /><br />2. On the same blog: "<a href="https://divisbyzero.com/2008/09/22/what-is-the-difference-between-a-theorem-a-lemma-and-a-corollary/">What is the difference between a theorem, a lemma, and a corollary?</a>"<br /><br />3. The "<a href="http://physicsinsights.org/glue_function.html">glue function</a>"<br /><br />4. Free <a href="http://linear.axler.net/LinearAbridged.html">abridged Linear Algebra</a> book from Sheldon Axler.</div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/caNPCNmgDn4" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/09/some-useful-math-links.htmltag:blogger.com,1999:blog-7379110960796014170.post-75871480611555653592017-09-16T11:06:00.000-04:002017-09-16T11:06:11.608-04:00Implicit Bias Test<div dir="ltr" style="text-align: left;" trbidi="on">I thoroughly enjoyed <a href="http://rationallyspeakingpodcast.org/show/rs-192-jesse-singal-on-the-problems-with-implicit-bias-tests.html">this</a> Jesse Singal interview on Rationally Speaking on the problems with the "<a href="https://en.wikipedia.org/wiki/Implicit-association_test">implicit association test</a>" for diagnosing implicit bias.<br /><br />The following <a href="https://www.youtube.com/watch?v=n5Q5FQfXZag">Dateline</a> video shows how the test was sold to the public as scientifically robust.<br /><br /><div class="separator" style="clear: both; text-align: center;"><iframe width="320" height="266" class="YOUTUBE-iframe-video" data-thumbnail-src="https://i.ytimg.com/vi/n5Q5FQfXZag/0.jpg" src="https://www.youtube.com/embed/n5Q5FQfXZag?feature=player_embedded" frameborder="0" allowfullscreen></iframe></div><br />For fun, you can <a href="https://implicit.harvard.edu/implicit/takeatest.html">take the test yourself</a>.<br /><br />For the problems with the test, check out Jesse Singal's piece from earlier this year, "<a href="http://nymag.com/scienceofus/2017/01/psychologys-racism-measuring-tool-isnt-up-to-the-job.html">Psychology’s Favorite Tool for Measuring Racism Isn’t Up to the Job</a>". It is a thoughtful essay, that should be read in its entirety.<blockquote class="tr_bq">A pile of scholarly work, some of it published in top psychology journals and most of it ignored by the media, suggests that the IAT falls far short of the quality-control standards normally expected of psychological instruments. The IAT, this research suggests, is a noisy, unreliable measure that correlates far too weakly with any real-world outcomes to be used to predict individuals’ behavior — even the test’s creators have now admitted as such. The history of the test suggests it was released to the public and excitedly publicized long before it had been fully validated in the rigorous, careful way normally demanded by the field of psychology.</blockquote>Singal is careful to point out that just because IAT is flawed it doesn't imply that implicit bias doesn't exist. I liked an analogy he used in the podcast. If a thermometer is flawed, you can't use it to determine if a person has a fever. The person may or may not have a fever, but the thermometer should probably be tossed away. </div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/gupZzfUJcWE" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/09/implicit-bias-test.htmltag:blogger.com,1999:blog-7379110960796014170.post-78944696895034369942017-09-12T12:39:00.000-04:002017-09-12T12:39:09.633-04:00Euler and Graph Theory<div dir="ltr" style="text-align: left;" trbidi="on">I have been enjoying Marcus du Sautoy's fine <a href="http://www.bbc.co.uk/programmes/b00srz5b/episodes/downloads">podcast</a> of famous mathematicians for BBC4. <div><br /></div><div>Yesterday, I listened to the Leonhard Euler episode. While I always knew Euler was one of the top mathematicians of all time, his <a href="https://en.wikipedia.org/wiki/Leonhard_Euler#Contributions_to_mathematics_and_physics">contributions</a> are truly remarkable.<div><br /></div><div>The podcast talks about how he solved the seven bridges of Konigsberg problem by inventing graph theory, and proving its first theorem. I looked at that theorem as it applies to a "kids game" in a <a href="http://sachinashanbhag.blogspot.com/2016/07/sunday-afternoon-fun.html">previous blog</a>.</div></div></div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/xKLgwBKkUIw" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/09/euler-and-graph-theory.htmltag:blogger.com,1999:blog-7379110960796014170.post-24559383829480851902017-09-05T18:31:00.000-04:002017-09-05T18:31:06.080-04:00Teacher's Day 2017<div dir="ltr" style="text-align: left;" trbidi="on">I did not expect writing this post would be so bittersweet. Last <a href="http://sachinashanbhag.blogspot.com/2016/09/teachers-day.html">Teacher's Day</a>, I decided I would use the occasion to highlight specific teachers, who have had an outsized impact on me.<br /><br />Today, I am going to tell you about Kartic C. Khilar, or KCK as he was called at IIT Bombay. KCK was a central figure, and participant, as I navigated a period of multiple transitions.<br /><br />Interestingly, I first "met" KCK even before I met him. The year I took the Joint Entrance Exam (JEE) to apply for admission to IIT, he was the principal administrator. The only reason I remember is because he had a "killer" last name (so juvenile, I know!).<br /><br />Like 200,000 other rats, I studied relentlessly for two years. JEE is like academic Olympics. We trained like mental athletes: cardio, weights, pilates, the whole nine yards. Then, the starting gun went off, and we scampered. The first two thousand got in.<br /><br />Miraculously, I tumbled my way into IIT Bombay first, and then to the chemical engineering department. KCK was the head of the department, when my "batch" arrived.<br /><br />He taught us fluid mechanics and solid-fluid operations. He was a fantastic teacher - one of the best I've had. His lectures were crisp. He was always cheerful. And he cared about all his students - not just toppers.<br /><br />He had one striking attribute: no ego. No made up sense of self-importance, which is all the more remarkable given the power gap between teachers and students (especially in India). If you went to his office, he would listen, despite how busy he was, or how unimportant you were.<br /><br />A highlight of the undergrad program at IIT is the B. Tech project (BTP), which is the undergrad equivalent of a PhD dissertation. Again, due to a random set of circumstances, he ended up being my BTP mentor. Over the course of the last year and half at IIT our interaction deepened, if only because we met one-on-one on a weekly basis to discuss research.<br /><br />Research in the Fluid Mechanics lab was fun. I don't think I would have embarked on a research career, if I hadn't enjoyed this experience so much. This work on "colloid-facilitated contaminant transport" with KCK and his grad student at that time - Tushar Sen - would end up becoming my first peer-reviewed <a href="http://www.sciencedirect.com/science/article/pii/S0927775703005545">publication</a>.<br /><br />I ended up at the University of Michigan as a grad student, in no small part due to his kind word. Michigan was his alma mater too. He visited Ann Arbor twice, while I was there. Once, when I was a PhD student, and later just before I started my new academic job at Florida State. Each time I went to Bombay, I would meet him; usually over lunch or dinner.<br /><br />Throughout this period, he selflessly offered his mind for me to pick, and his ocean of experience for me to draw from. At several points during this journey, I abandoned hopes of an academic career. Each time, he listened without judgment, and quietly held a mirror to my desire for autonomy and passion for teaching. For better or for worse, he was instrumental in me ending up on the trajectory I am currently on.<br /><br />And I couldn't be more grateful! Sometimes you try to peek over the horizon, but you can't see what a taller person who has been to more places can (in my case, that is literally true too).<br /><br /><br />In 2009, I shut the door to my office and wept, when I learnt about his <a href="http://www.che.iitb.ac.in/online/story/obituary-prof-kartic-khilar-a-tribute-his-career-and-compassion">untimely passing</a>. He was 57, in great mental and physical shape, and I always expected him to be around forever.<br /><br />When I first encountered KCK in 1994, I knew him as an administrator. Later at IIT he became my chairman and teacher, before becoming my BTP supervisor.<br /><br />Somewhere along the way, he became a mentor, and a close friend; emails that started with "Dear Prof. Khilar" eventually started with "Dear Kartic".<br /><br />Today, even though I knew it would bounce, I nearly wrote (to his familiar email address), "Dear Kartic, you are sorely missed." </div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/gk3DA0q1xvY" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/09/teachers-day-2017.htmltag:blogger.com,1999:blog-7379110960796014170.post-84800244643808376002017-09-02T18:29:00.002-04:002017-09-02T18:29:50.342-04:00Complex Numbers: Part Deux<div dir="ltr" style="text-align: left;" trbidi="on">I was pointed to this excellent <a href="https://www.youtube.com/playlist?list=PLiaHhY2iBX9g6KIvZ_703G3KJXapKkNaF">series on complex numbers from Welch labs</a>, following my last post on complex numbers. It in the <a href="https://www.youtube.com/channel/UCYO_jab_esuFRV4b17AJtAw">3Blue1Brown</a> mold, with just the right dose of insight and animation. The complex number series starts with basic ideas, and ends with a discussion of Riemann surfaces.<br /><br />I also came across an interesting way of proving exp(ix) = cos x + i sin x (@fermatslibrary), which I feel compelled to share, since we are already talking about complex numbers.<br /><br />Let \(f(x) = e^{-ix} (\cos x + i \sin x)\).<br /><br />The derivative of this function is \[f'(x) = e^{-ix} (i\cos x - i \sin x) - i e^{-ix} (\cos x + i \sin x) = 0.\] Since \(f'(x) = 0\), the function is a constant.<br /><br />Also f(0) = 1, which implies f(x) = 1.<br /><br />Thus, \(e^{ix} = \cos x + i \sin x\).<br /><br /><b>PS</b>: One of my students told me last week about the new podcast (<a href="https://www.benbenandblue.com/">Ben, Ben, and Blue</a>) that Grant Sanderson (of 3Blue1Brown) hosts on math, computer science and education. It is delightful.</div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/b_lUtdfBtg0" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/09/complex-numbers-part-deux.htmltag:blogger.com,1999:blog-7379110960796014170.post-80034317707909862572017-08-29T20:27:00.000-04:002017-08-29T20:27:01.433-04:00Anomalous Diffusion<div dir="ltr" style="text-align: left;" trbidi="on">I've been taking a deep dive into the world of anomalous diffusion over the past month. It is a fascinating subject that integrates applications from a variety of different fields.<br /><br />For someone interested, I'd recommend the following resources:<br /><br />1. A Physics World feature on "Anomalous diffusion spreads its wings" (<a href="http://www.tau.ac.il/~klafter1/ar1.pdf">pdf</a> - currently not paywalled)<br /><br />2. A <a href="https://www.youtube.com/watch?v=ZKjQKWq02_4">YouTube video</a> on anomalous diffusion in crowded environments<br /><br /><div class="separator" style="clear: both; text-align: center;"><iframe allowfullscreen="" class="YOUTUBE-iframe-video" data-thumbnail-src="https://i.ytimg.com/vi/ZKjQKWq02_4/0.jpg" frameborder="0" height="266" src="https://www.youtube.com/embed/ZKjQKWq02_4?feature=player_embedded" width="320"></iframe></div><br /><br />3. A gentle introduction/tutorial on <a href="https://arxiv.org/pdf/0805.0419.pdf">normal and anomalous diffusion</a>, which introduces the intuition and mechanics of fractional calculus<br /><br />4. A more academic <a href="http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.473.7695&rep=rep1&type=pdf">review</a> of anomalous diffusion and fractional dynamics (may be <a href="http://www.sciencedirect.com/science/article/pii/S0370157300000703">paywalled</a>)</div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/HWzmQrddno0" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/08/anomalous-diffusion.htmltag:blogger.com,1999:blog-7379110960796014170.post-22103573058764484342017-08-23T10:10:00.003-04:002017-08-23T10:10:50.136-04:00If $1 = 1 sec ...<div dir="ltr" style="text-align: left;" trbidi="on">If $1 were equal to 1 second, the median US household income per year of $50,000 would correspond to half a day.<br /><br />This helps puts millions, billions, and trillions into perspective.<br /><br />Roughly,<br /><ul style="text-align: left;"><li>$1 million = 2 weeks</li><li>$1 billion = 32 years</li><li>$1 trillion = 300 centuries (before recorded history)</li></ul>A trillion is a really large number! </div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/IHZqc2SN44E" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/08/if-1-1-sec.htmltag:blogger.com,1999:blog-7379110960796014170.post-69677843180920213992017-08-15T12:49:00.002-04:002017-08-15T12:49:40.166-04:00Diffusion: A Historical Perspective<div dir="ltr" style="text-align: left;" trbidi="on">The paper (<a href="https://pdfs.semanticscholar.org/24c0/5e77599f7ec211b9c0bbf326138607889415.pdf">pdf</a>) "One and a half century of diffusion: Fick, Einstein, Before and Beyond" by Jean Philibert traces the history of diffusion phenomena.<br /><br />It starts with <a href="https://en.wikipedia.org/wiki/Thomas_Graham_(chemist)">Thomas Graham</a> (of dialysis fame) who perhaps made the first systematic observations, which were integrated into phenomenological law by German physiologist <a href="https://en.wikipedia.org/wiki/Adolf_Eugen_Fick">Adolf Fick</a> in 1855, at the age of 26.<br /><br />Fick observed the analogy between mass diffusion and heat conduction (now considered obvious), and piggy-backed on Fourier's law of conduction (1822). The paper cites the opening lines of Fick's work:<br /><blockquote class="tr_bq">A few years ago, Graham published an extensive investigation on the diffusion of salts in water, in which he more especially compared the diffusibility of different salts. It appears to me a matter of regret, however, that in such an exceedingly valuable and extensive investigation, the development of a fundamental law, for the operation of diffusion in a single element of space, was neglected, and I have therefore endeavoured to supply this omission.</blockquote>Next, the paper talks about the contributions of W. C. Roberts-Austen (an assistant to Thomas Graham, and successor as Master of the Mint) to quantification of diffusion in solids.<br /><br />In 1905, Einstein integrated Robert Brown's observations of random zig-zag trajectories and Fick's phenomenological laws, with the crucial observation that it was the mean-squared displacement, and not the mean displacement that was related to diffusion.<br /><br />Following Einstein's paper, the experimental work of Perrin was responsible helping the world accept the link between the microscopic (MSD is proportional to diffusivity and time) and macroscopic worlds (flux is proportional to concentration gradient).<br /><br />It is always interesting to look at the chronological development of (now familiar) ideas. These uncontroversial ideas were once strongly wrestled with. It took centuries for scientists to come up with a comprehensive understanding, and to develop interesting applications based off of it.</div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/RlwSUOuAfIg" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com2http://sachinashanbhag.blogspot.com/2017/08/diffusion-historical-perspective.htmltag:blogger.com,1999:blog-7379110960796014170.post-67709878068357351832017-08-12T14:30:00.000-04:002017-08-12T14:31:00.383-04:00Exam Question on Fitting Sums of Exponentials to Data<div dir="ltr" style="text-align: left;" trbidi="on">I wrote the question below for our PhD qualifiers. It addresses a problem I have been thinking about for over a decade now - starting from my time as a graduate student: how to fit a sum of decaying exponentials to data?<br /><br />The question explores a method called the <a href="https://en.wikipedia.org/wiki/Prony%27s_method">Prony</a> method. Here is the question:<br /><br />A classical problem in data analysis involves fitting a sum of exponentials to a time series of uniformly sampled observations. Here, let us suppose we are given N observations \((t_i, f_i)\), where \(t_i = i \Delta t\) for \(i = 0, 1, ..., N-1\).<br /><br />We want to fit the data to a sum of two exponentials. The <b>model equation</b> is, \[\hat{f}(t) = a_1 e^{b_1 t} + a_2 e^{b_2 t}.\] The general nonlinear regression problem to determine \(\{a_j, b_j\}\) becomes difficult as the number of exponentials in the sum increases. A number of quasi-linear methods have been developed to address this. In the question, we will explore one of these methods, and determine the fitting parameters.<br /><div><br />(a) First, generate a synthetic dataset \((t_i, f_i)\) with true \(a_1^* = a_2^* = 1.0\), \(b_1^* = -2.0\), \(b_2^* = -0.2\). Use \(t_0 = 0\), \(\Delta t = 1\), and N = 20. Attach a plot of the synthetic dataset. Use this dataset for numerical calculations below.<br /><br />(b) If \(b_1\) and \(b_2\) are known, then we can determine \(a_1\) and \(a_2\) by linear least squares. Set \(u_1 = e^{b_1 \Delta t}\) and \(u_2 = e^{b_2 \Delta t}\). Recognize that \(e^{b_i t_j} = e^{b_i j \Delta t} = u_i^j\). Hence from the model eqn, we can get a <b>linear system</b>:<br />\begin{align}<br />f_0 & = a_1 u_1^0 + a_2 u_2^0 \nonumber\\<br />f_1 & = a_1 u_1^1 + a_2 u_2^1 \nonumber\\<br />\vdots & = \vdots \nonumber\\<br />f_{N-1} & = a_1 u_1^{N-1} + a_2 u_2^{N-1}<br />\end{align}<br />Write a program to determine \(a_1\) and \(a_2\), given the data, \(b_1\) and \(b_2\).<br /><br />(c) Consider the polynomial \(p(z)\), which has \(u_1\) and \(u_2\) as its roots, \(p(z) = (z-u_1)(z-u_2) = z^2 - d_1 z -d_2 = 0\). Express \(u_1\) and \(u_2\) in terms of \(d_1\) and \(d_2\).<br /><br />(d) Now we seek to take linear combinations equations in the linear system above with the goal of eliminating \(a_j\). For example, consider the first three equations. If we multiply the first eqn by \(d_2\), the next by \(d_1\), and the third by -1 and sum them up.<br />\begin{align*}<br />d_2 f_0 & = a_1 d_2 + a_2 d_2\\<br />d_1 f_1 & = a_1 u_1 d_1 + a_2 u_2 d_1 \\<br />-1 f_2 & = -a_1 u_1^2 - a_2 u_2^2.<br />\end{align*}<br />We get \(-F_2 +d_1 F_1 + d_2 F_0 = -a_1(u_1^2 - d_1 u_1 - d_2) -\) \( a_2(u_2^2 -d_1 u_2 - d_2) = 0\), since \(p(u_i) = 0\).</div><div><br />We can pick the next set of three equations, and repeat the process (multiply by \(d_2\), \(d_1\), and -1 before summing up). Show that we end up with the following linear system:<br />\[\begin{bmatrix} f_{1} & f_0 \\ f_2 & f_1 \\ <br />\vdots & \vdots \\<br />f_{N-2} & f_{N-3} \\<br />\end{bmatrix} \begin{bmatrix} d_1 \\ d_2 \end{bmatrix} = \begin{bmatrix} f_2 \\ f_{3} \\ \vdots \\ f_{N-1} \end{bmatrix}\]<br />Determine \(d_1\) and \(d_2\), and hence \(u_1\) and \(u_2\). From this, find the estimated \(b_1\) and \(b_2\).<br /><br />(e) Once you know \(b_1\) and \(b_2\) find \(a_1\) and \(a_2\) by linear least squares solution of linear system.<br /><br /></div></div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/tTNGdIInFrM" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/08/exam-question-on-fitting-sums-of.htmltag:blogger.com,1999:blog-7379110960796014170.post-10001900818720983522017-08-09T16:35:00.000-04:002017-08-09T23:05:02.175-04:00Complex Numbers<div dir="ltr" style="text-align: left;" trbidi="on">1. What really are <a href="http://robjlow.blogspot.com/2017/06/what-is-this-thing-called-i.html?m=1">complex numbers</a>?<br /><br />2. The joy of <a href="https://mathwithbaddrawings.com/2017/05/17/the-joy-of-slightly-fishy-proofs/">slightly fishy proofs</a>.<br /><br />3. This <a href="https://math.stackexchange.com/questions/4961/interesting-results-easily-achieved-using-complex-numbers">discussion</a> on MathOverflow</div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/DWrFECnrRo4" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com3http://sachinashanbhag.blogspot.com/2017/08/complex-numbers.htmltag:blogger.com,1999:blog-7379110960796014170.post-24432746365791598692017-08-02T15:52:00.000-04:002017-08-02T15:52:04.594-04:00NumPy and Matlab<div dir="ltr" style="text-align: left;" trbidi="on">This post bookmarks two sites that provide handy cheat sheets of numpy equivalents for Matlab/Octave commands.<br /><br />The ones for linear algebra are particularly handy, because that is one subdomain where Matlab's notation is more natural.<br /><br />1. Numpy for Matlab users (<a href="http://mathesaurus.sourceforge.net/matlab-numpy.html">Mathesaurus</a>)<br /><br />2. Cheatsheets for Numpy, Matlab, and Julia (<a href="https://cheatsheets.quantecon.org/">quantecon</a>)</div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/vcft0PJS-hs" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/08/numpy-and-matlab.htmltag:blogger.com,1999:blog-7379110960796014170.post-49677357244099734842017-07-28T21:19:00.000-04:002017-07-24T08:59:15.425-04:00Interesting Scaling Laws<div dir="ltr" style="text-align: left;" trbidi="on">I recently read Geoffrey West's book "<a href="https://www.amazon.com/Scale-Universal-Innovation-Sustainability-Organisms/dp/1594205582">Scale</a>", and thought it was really great. Here are some resources to prime you for the subject.<br /><br />1. <a href="https://www.ted.com/talks/geoffrey_west_the_surprising_math_of_cities_and_corporations">TED Talk</a><br /><br />2. Talk <a href="https://www.youtube.com/watch?v=GoHD1ROPiUc">@ Google</a><br /><br />3. Essay at the <a href="https://www.edge.org/conversation/geoffrey_west-why-cities-keep-growing-corporations-and-people-always-die-and-life-gets">Edge</a><br /><br />4. Essay on <a href="https://medium.com/sfi-30-foundations-frontiers/scaling-the-surprising-mathematics-of-life-and-civilization-49ee18640a8">Medium</a></div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/vZqgD0Jj11I" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/06/interesting-scaling-laws.htmltag:blogger.com,1999:blog-7379110960796014170.post-41031123725894969792017-07-25T10:44:00.000-04:002017-07-25T10:44:00.167-04:00Russell's paradox<div dir="ltr" style="text-align: left;" trbidi="on"><div class="tr_bq">I came across this interesting paradox on a recent podcast. According to <a href="https://en.wikipedia.org/wiki/Russell%27s_paradox">wikipedia</a>:</div><blockquote>According to naive set theory, any definable collection is a set. Let ''R'' be the set of <b>all</b> sets that are not members of themselves. If ''R'' is not a member of itself, then its definition dictates that it must contain itself, and if it contains itself, then it contradicts its own definition as the set of all sets that are not members of themselves. This contradiction is Russell's paradox. </blockquote><blockquote>Symbolically:<br />\[\text{Let } R = \{ x \mid x \not \in x \} \text{, then } R \in R \iff R \not \in R\]</blockquote>There is a nice commentary on the paradox in <a href="https://www.scientificamerican.com/article/what-is-russells-paradox/">SciAm</a>, and a superb entry on the <a href="https://plato.stanford.edu/entries/russell-paradox/">Stanford Encyclopedia of Philosophy</a></div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/y5s1tE1Vc8o" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/07/russells-paradox.htmltag:blogger.com,1999:blog-7379110960796014170.post-38426202088200803172017-07-19T10:48:00.000-04:002017-07-19T10:48:20.013-04:00Questions Kids Ask<div dir="ltr" style="text-align: left;" trbidi="on">Between my curious 4- and 8-year olds, I got asked the following questions in the past month.<br /><br />I found all of them fascinating.<br /><br />1. Why are our front milk teeth (incisors) the first to fall out?<br />2. Why is "infinity minus infinity" not equal to zero?<br />3. Why don't you get a rainbow when you shine a flashlight on rain in the night?<br />4. How are Cheerios and donuts made (into tori)?<br />5. His, hers, ours, yours. Then why not "mines"?<br /><br />PS: I also learned from my 4-year old that <a href="http://spiders.ucr.edu/daddylonglegs.html">daddy long legs aren't really spiders</a> and don't spin webs, and that <a href="https://www.google.com/search?q=sea+turtles+jellyfish&source=lnms&tbm=isch&sa=X&ved=0ahUKEwiPobTbg47VAhVEOSYKHckgB9kQ_AUICigB&biw=1920&bih=963">sea turtles feed on jellyfish</a>.</div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/ecIMR2JPj8g" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/07/questions-kids-ask.htmltag:blogger.com,1999:blog-7379110960796014170.post-73248781229028323952017-07-16T10:37:00.000-04:002017-07-16T10:37:02.284-04:00Matplotlib: Subplots, Inset Plots, and Twin Y-axes<div dir="ltr" style="text-align: left;" trbidi="on">This <a href="https://gist.github.com/shane5ul/9e8dbade9f9f9f4de5c8b00a41b53f20">jupyter notebook</a> highlights ways in which matplotlib gives you control over the layout of your charts. This is intended as a personal cheatsheet.<br /><div><br /></div><script src="https://gist.github.com/shane5ul/9e8dbade9f9f9f4de5c8b00a41b53f20.js"></script> <br /><div><br /></div></div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/ibswlrsV8fY" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com1http://sachinashanbhag.blogspot.com/2017/07/matplotlib-subplots-inset-plots-and.htmltag:blogger.com,1999:blog-7379110960796014170.post-84514712745093005092017-07-07T16:34:00.002-04:002017-07-07T16:34:16.832-04:00John Roberts Commencement Speech<div dir="ltr" style="text-align: left;" trbidi="on"><div class="separator" style="clear: both; text-align: left;">This part of the address is really nice and timeless.</div><div class="separator" style="clear: both; text-align: left;"><br /></div><div class="separator" style="clear: both; text-align: center;"><iframe width="320" height="266" class="YOUTUBE-iframe-video" data-thumbnail-src="https://i.ytimg.com/vi/pavcM8VGtX0/0.jpg" src="https://www.youtube.com/embed/pavcM8VGtX0?feature=player_embedded" frameborder="0" allowfullscreen></iframe></div><br />The transcript of the full speech is available <a href="http://time.com/4845150/chief-justice-john-roberts-commencement-speech-transcript/">here</a>.</div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/7spx4ii8BcI" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/07/john-roberts-commencement-speech.htmltag:blogger.com,1999:blog-7379110960796014170.post-84717300133956368502017-07-05T16:20:00.001-04:002017-07-05T16:20:07.674-04:00Joints from Marginals: Compilation<div dir="ltr" style="text-align: left;" trbidi="on">For convenience, here is a link to the three blogs in this series in one place.<div><br /></div><div>1. A <a href="https://sachinashanbhag.blogspot.com/2017/06/joint-distribution-from-marginals.html">technique</a> for solving the problem in a special case</div><div><br /></div><div>2. The <a href="https://sachinashanbhag.blogspot.com/2017/06/joint-from-marginals-why.html">reason</a> this technique works</div><div><br /></div><div>3. The <a href="https://sachinashanbhag.blogspot.com/2017/07/joint-from-marginals-non-gaussian.html">corners/edges</a> of this technique, or how it fails for non-Gaussian marginals</div></div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/kek9wb_HrT4" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/07/joints-from-marginals-compilation.htmltag:blogger.com,1999:blog-7379110960796014170.post-50122949468697229432017-07-02T09:32:00.001-04:002017-07-02T09:32:09.644-04:00Joint from Marginals: non-Gaussian Marginals<div dir="ltr" style="text-align: left;" trbidi="on">In a <a href="http://sachinashanbhag.blogspot.com/2017/06/joint-from-marginals-why.html">previous</a> post, I asked the question if the method described here can be used with non-Gaussian distributions.<br /><br />Let us explore that by considering two independent zero mean, unit variance distributions that are not Gaussian. Let us sample \(x_1\) from a <a href="https://en.wikipedia.org/wiki/Triangular_distribution">triangular distribution</a>, and \(x_2\) from a <a href="https://en.wikipedia.org/wiki/Uniform_distribution_(continuous)">uniform distribution</a>.<br /><br />We consider a triangular distribution with zero mean and unit variance, which is symmetric about zero (spans -sqrt(6) to +sqrt(6)). Similarly, we consider a symmetric uniform distribution, which spans -sqrt(3) to +sqrt(3).<br /><br />Samples from these independent random variables are shown below.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://4.bp.blogspot.com/-pNLDuhtXQ5Q/WRyrwnETwYI/AAAAAAAADCM/pUpMWJIS4Lk-lt8_NU8mW6-v-H2pPpd3gCLcB/s1600/1.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="400" src="https://4.bp.blogspot.com/-pNLDuhtXQ5Q/WRyrwnETwYI/AAAAAAAADCM/pUpMWJIS4Lk-lt8_NU8mW6-v-H2pPpd3gCLcB/s400/1.png" width="400" /></a></div><div class="separator" style="clear: both; text-align: left;">When we use a correlation coefficient of 0.2, and use the previous recipe, we get correlated random variables with zero mean and the same covariance matrix, but ...</div><div class="separator" style="clear: both; text-align: center;"><a href="https://3.bp.blogspot.com/--bINPDVldFQ/WRysK_14TvI/AAAAAAAADCQ/eaKLJTfyNxkH5_QBSP3sfULOUPgnsa91gCLcB/s1600/2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="400" src="https://3.bp.blogspot.com/--bINPDVldFQ/WRysK_14TvI/AAAAAAAADCQ/eaKLJTfyNxkH5_QBSP3sfULOUPgnsa91gCLcB/s400/2.png" width="400" /></a></div><div class="separator" style="clear: both; text-align: left;">... the marginals are not exactly the same!</div><div class="separator" style="clear: both; text-align: left;"><br /></div><div class="separator" style="clear: both; text-align: left;">This is evident when we increase the correlation coefficient to say 0.5.</div><div class="separator" style="clear: both; text-align: left;"><br /></div><div class="separator" style="clear: both; text-align: center;"><a href="https://2.bp.blogspot.com/-4dm9JM45UO4/WRysg793OmI/AAAAAAAADCU/8jZMoOYiNs0jLXr8g7nQYdpAC5lZYnRzACLcB/s1600/2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="400" src="https://2.bp.blogspot.com/-4dm9JM45UO4/WRysg793OmI/AAAAAAAADCU/8jZMoOYiNs0jLXr8g7nQYdpAC5lZYnRzACLcB/s400/2.png" width="400" /></a></div><div class="separator" style="clear: both; text-align: left;">The sharp edges of the uniform distribution get smoothened out.</div><div class="separator" style="clear: both; text-align: left;"><br /></div><div class="separator" style="clear: both; text-align: left;">Did the method fail?</div><div class="separator" style="clear: both; text-align: left;"><br /></div><div class="separator" style="clear: both; text-align: left;">Not really. If you paid attention, the method is designed to preserve the mean and the covariance matrix (which is does). It doesn't really guarantee the preservation of the marginal distributions. </div><br /></div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/ya05YBWZAU0" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/07/joint-from-marginals-non-gaussian.htmltag:blogger.com,1999:blog-7379110960796014170.post-7025669095900254742017-06-28T12:38:00.000-04:002017-06-28T12:38:33.872-04:00Printing webpages as PDFs<div dir="ltr" style="text-align: left;" trbidi="on"><a href="https://www.printfriendly.com/">PrintFriendly and PDF</a> has a useful browser extension (tested on Chrome) that creates more readable PDFs from web content.<br /><br />Here is a screenshot (click to enlarge) from a Matlab blog that I follow:<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://4.bp.blogspot.com/-JKuxK6ID6fg/WRXjeSwHWaI/AAAAAAAADA0/URIxezSfYH88BbvoEDI-8ikSmhXCJF9MwCEw/s1600/Screenshot%2Bfrom%2B2017-05-12%2B12-25-09.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="320" src="https://4.bp.blogspot.com/-JKuxK6ID6fg/WRXjeSwHWaI/AAAAAAAADA0/URIxezSfYH88BbvoEDI-8ikSmhXCJF9MwCEw/s400/Screenshot%2Bfrom%2B2017-05-12%2B12-25-09.png" width="400" /></a></div><br /><br /><div class="separator" style="clear: both; text-align: center;"></div>Notice that the webpage has lots of links, and a frame on the left.<br /><br />When I use the "Print to File" feature directly from my Chrome browser, I get a PDF which looks like this:<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-6cvqkJ4aseg/WRXjhC5uYvI/AAAAAAAADBA/wvqXVdLx08kmeeJYCdGCC4_8PPUGCYMvwCEw/s1600/Screenshot%2Bfrom%2B2017-05-12%2B12-25-48.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="347" src="https://1.bp.blogspot.com/-6cvqkJ4aseg/WRXjhC5uYvI/AAAAAAAADBA/wvqXVdLx08kmeeJYCdGCC4_8PPUGCYMvwCEw/s400/Screenshot%2Bfrom%2B2017-05-12%2B12-25-48.png" width="400" /></a></div>It does the job, but it looks very amateurish. On more complicated websites, results can be horrendous.<br /><br />Here is the same webpage, now using PrintFriendly.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://4.bp.blogspot.com/-CfnI9BC8rMk/WRXjrwdt01I/AAAAAAAADBA/7zEj50gbvO4KwYnK3ukdeSyHpY3pc5qLgCEw/s1600/Screenshot%2Bfrom%2B2017-05-12%2B12-26-57.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="400" src="https://4.bp.blogspot.com/-CfnI9BC8rMk/WRXjrwdt01I/AAAAAAAADBA/7zEj50gbvO4KwYnK3ukdeSyHpY3pc5qLgCEw/s400/Screenshot%2Bfrom%2B2017-05-12%2B12-26-57.png" width="307" /></a></div>Notice that the PDF is much cleaner, is well formatted, and contains all the relevant information.</div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/xK-TqnptFVs" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/06/printing-webpages-as-pdfs.htmltag:blogger.com,1999:blog-7379110960796014170.post-83125201019034120102017-06-22T15:53:00.000-04:002017-06-22T15:53:20.650-04:00Joint from Marginals: Why?<div dir="ltr" style="text-align: left;" trbidi="on">In the <a href="http://sachinashanbhag.blogspot.com/2017/06/joint-distribution-from-marginals.html">previous</a> blog post, we saw a special example in which we were able to sample random variables from a joint 2D-Gaussian distribution from the marginals and the correlation coefficient.<br /><br />I listed a simple method, which seemed to work like magic. It had two simple steps:<br /><br /><ul style="text-align: left;"><li>Cholesky decomposition of the covariance matrix, C(Y)</li><li>Y = LX, where X are independent random variables</li></ul><br />The question is, why did the method work?<br /><br />Note that the covariance matrix of random variables with zero mean and unit standard deviation can be written as, \(C(Y) = E(Y Y')\), where \(E()\) denotes the expected value of a random variable. Thus, we can write the expected value of the Y generated by the method as, \[\begin{align*} E(Y Y') & = E\left(LX (LX)'\right)\\ & = L E(XX') L' \\ & = L I L'\\ & = LL' = C.\end{align*}.\] Here we used the fact that the covariance of X is an identity matrix by design.<br /><br />Note that this method preserves the covariance matrix (and hence the standard deviation of the marginals).<br /><br />Does it preserve the mean?<br /><br />Yes. \(E(Y) = E(LX) = L E(X) = 0.\)<br /><br />Do the marginals have to be normal for this method to work? Would this work for any distribution (with zero mean, and unit standard deviation)?<br /><br />We will explore this in a subsequent blog.<br /><br /></div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/jIIUr8RAKio" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com0http://sachinashanbhag.blogspot.com/2017/06/joint-from-marginals-why.htmltag:blogger.com,1999:blog-7379110960796014170.post-75074865068231560232017-06-15T12:06:00.001-04:002017-06-15T12:06:48.930-04:00Joint Distribution From Marginals<div dir="ltr" style="text-align: left;" trbidi="on">Consider two dependent random variables, \(y_1\) and \(y_2\), with a correlation coefficient \(\rho\).<br /><br />Suppose you are given the marginal distributions \(\pi(y_1)\) and \(\pi(y_2)\) of the two random variables. Is it possible to construct the joint probability distribution \(\pi(y_1, y_2)\) from the marginals?<br /><div><br /></div>In general, the answer is no. There is no unique answer. The marginals are like shadows of a hill from two orthogonal angles. The shadows are not sufficient to specify the full 3D shape (joint distribution) of the hill.<br /><br />Let us simplify the problem a little, so that we can seek a solution.<br /><br />Let us assume \(y_1\) and \(y_2\) have zero mean and unit standard deviation. We can always generalize later by shifting (different mean) and scaling (different standard distribution). Let us also stack them into a single random vector \(Y = [y_1, y_2]\).<br /><br />The <a href="https://en.wikipedia.org/wiki/Covariance_matrix">covariance matrix</a> of two such random variables is given by, \[C(Y) = \begin{bmatrix} E(y_1 y_1) - \mu_1 \mu_1 & E(y_1 y_2) - \mu_1 \mu_2 \\ E(y_2 y_1) - \mu_2 \mu_1 & E(y_2 y_2) - \mu_2 \mu_2 \end{bmatrix} = \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix},\] where \(\mu\) and \(\sigma\) refer to the mean and standard deviation.<br /><div style="text-align: left;"><br /></div><div style="text-align: left;"><b>Method</b></div><br />A particular method for sampling from the joint distribution of correlated random variables \(Y\) begins by drawing samples of independent random variables \(X = [x_1, x_2]\) which have the same distribution as the desired marginal distributions.<br /><br />Note that the covariance matrix in this case is an identity matrix, because the correlation between independent variables is zero \(C(X) = I\).<br /><br />Now we recognize that the covariance matrix \(C(Y)\) is symmetric and positive definite. We can use Cholesky decomposition \(C(Y) = LL^T\) to find the lower triangular matrix \(L\).<br /><br />The recipe then says that we can draw the correlated random variables with the desired marginal distribution by simply setting \(Y = L X\).<br /><br /><b>Example</b><br /><b><br /></b>Suppose we seek two random variables whose marginals are normal distributions (zero mean, unit standard deviation) with a correlation coefficient 0.2.<br /><br />The method above asks us to start with independent random variables \(X\) such as those below.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://4.bp.blogspot.com/-yc6SE7nazqM/WRylDZ7nDOI/AAAAAAAADBw/bOsSEhSFXjgPTN0N5AgypQ0xWBO4aceqwCLcB/s1600/1.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="400" src="https://4.bp.blogspot.com/-yc6SE7nazqM/WRylDZ7nDOI/AAAAAAAADBw/bOsSEhSFXjgPTN0N5AgypQ0xWBO4aceqwCLcB/s400/1.png" width="400" /></a></div>Cholesky decomposition with \(\rho\) = 0.2, gives us, \[L = \begin{bmatrix} 1 & 0 \\ 0.1 & 0.9797 \end{bmatrix}.\] If we generate \(Y = LX\) using the same data-points used to create the scatterplot above, we get,<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://3.bp.blogspot.com/-MB0nLXdDD9s/WRymodTmzmI/AAAAAAAADB8/lEaz4NaCekE9h8AQ7iL8BtNHLzgfZeSdACLcB/s1600/2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="400" src="https://3.bp.blogspot.com/-MB0nLXdDD9s/WRymodTmzmI/AAAAAAAADB8/lEaz4NaCekE9h8AQ7iL8BtNHLzgfZeSdACLcB/s400/2.png" width="400" /></a></div>It has the same marginal distribution, and a non-zero correlation coefficient as is visible from the figure above.</div><img src="http://feeds.feedburner.com/~r/CluelessFundatma/~4/qPNg-uk_RDM" height="1" width="1" alt=""/>Sachin Shanbhaghttps://plus.google.com/115150474038005608083noreply@blogger.com1http://sachinashanbhag.blogspot.com/2017/06/joint-distribution-from-marginals.html