tag:blogger.com,1999:blog-44538050959428128632021-07-16T21:16:27.303-07:00Surprise and Coincidence - musings from the long tailTed Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.comBlogger55125tag:blogger.com,1999:blog-4453805095942812863.post-44230672813696004792021-06-07T09:03:00.007-07:002021-06-07T16:41:27.265-07:00What's new in T-Digest New Version 3.3?<p>There's a new version of $t$-digest in town. This new version doesn't much change how things work from the standpoint of the API, but it does improve accuracy under most circumstances. If you are using $t$-digest, it is a generally good idea to upgrade to the 3.2 release.</p><p>Let's dig into the changes after a quick review of how the $t$-digest works.</p><h3 style="text-align: left;">Density estimation through clustering</h3><div>The basic idea behind the $t$-digest is that if you cluster your data, the position and number of samples in each cluster can be used to estimate the distribution of your original data. Further, if you guarantee that the clusters near the extreme values are kept small and the clusters near the median are allowed to grow, you can maintain very good accuracy for extreme quantiles while sacrificing only a small amount of accuracy for central values. Because clustering in one dimension is very different from clustering in multiple dimensions, the $t$-digest can provide a very good trade-off between memory consumption and accuracy.</div><div><br /></div><div>The $t$-digest guarantees not to use more than a set amount of memory, but this can theoretically be shown to make it impossible to give guarantees on the relative accuracy of quantile estimates. Indeed, a pathological distribution of data can cause a $t$-digest to give large errors. Happily, a sufficiently pathological distribution of values is essentially impossible to have in any physical measurement. For instance, if you have times that range from picoseconds to the age of the universe with massive skew, you will still be hundreds of powers of ten short of a distribution sufficiently pathological to cause large errors.</div><div><br /></div><div>The sizes of the clusters in a $t$-digest are bounded using a so-called scale function and a compression factor. The scale factor determines how aggressive the digest is about keeping extremal clusters small and the compression factor limits the number of clusters retained. </div><h3 style="text-align: left;">Better boundary conditions</h3><p>The latest version of the $t$-digest maintains several invariants more assiduously than previous version. First of all, the minimum and maximum values ever seen are retained explicitly. Secondly, the digest makes sure that the first and last centroid always have unit weight except when the K_0 scale function is used. This improves the quality of interpolation at the end points as we will see in the next section.</p><h3 style="text-align: left;">Better interpolation</h3><p>The most important accuracy changes to the $t$-digest (besides general code simplification and cleanups) have to do with the way that interpolation is handled. The basic idea is that centroids that represent a single point need special handling because we know exactly where that point is. Moreover, we know for the first and last centroid that at least one sample is at the global minimum (or maximum respectively) and thus we can treat the first (or last) centroid as if it has one fewer sample than observed. With the K_0 scale function, extreme centroids can have more than one sample and if such an extreme centroid has exactly two samples, we can infer the exact location of the second sample allowing us to do better interpolation. Even if there are more samples than that, we can pull the extreme sample out of the mean and get slightly better interpolation.</p><p>The case where a singleton centroid is next to a centroid with more samples is handled better as well.</p><p>How this works is illustrated here. Here we have two cases with three centroids. On the left, all three centroids have weight greater than one while on the right, two centroids have weight equal to one.</p><p></p><div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-jP8CgSPF7xI/YLwdVlGkZII/AAAAAAAAyi0/Y04nybQJZ9c4RcX7KSJndwZ_sjdYwZulwCLcBGAsYHQ/s2163/combined.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="663" data-original-width="2163" height="137" src="https://1.bp.blogspot.com/-jP8CgSPF7xI/YLwdVlGkZII/AAAAAAAAyi0/Y04nybQJZ9c4RcX7KSJndwZ_sjdYwZulwCLcBGAsYHQ/w445-h137/combined.png" width="445" /></a></div><p style="text-align: left;">The basic idea is that we normally assume that half of the points for each centroid land on either side of the centroid. That means that the number of samples between two centroids is assumed, for interpolation purposes, to be the sum of the halves from the two adjacent centroids. But when we have a neighboring centroid with just one sample, we don't have to linearly interpolate all the way up to that centroid and thus can spread the weight from the non-singleton centroid over a slightly narrower range.</p><p>This is important near the tails of any digest, but it is really important to get this right when the total amount of data that has been added to a digest is a small-ish multiple of the compression factor of the digest. This is because any such digest will have many singleton centroids near the tail and a clean transition to non-singleton centroids near the center can be a major source of error.</p><p>The impact of this can be seen in comparisons against the relative error version of the KLL digest. Here the compression factor of the $t$-digest was adjusted to give a stored digest size of roughly 1kB or 4kB. The first value represents relatively normal operation for a $t$-digest with compression factor in the range of 100 or so while the larger value matches the size of the KLL digest. The following figure shows this comparison for varying numbers of samples.</p><p></p><div class="separator" style="clear: both; text-align: center;"><div class="separator" style="clear: both; text-align: center;"><a href="https://lh3.googleusercontent.com/-XeWnKFhNMXI/YLwgIqZOG-I/AAAAAAAAyjA/NvQKXc3giFImV3Bb3d6FGo3IYiP3KXztwCLcBGAsYHQ/kll-comparison-x.png" style="margin-left: 1em; margin-right: 1em;"><img alt="" data-original-height="1227" data-original-width="2048" height="281" src="https://lh3.googleusercontent.com/-XeWnKFhNMXI/YLwgIqZOG-I/AAAAAAAAyjA/NvQKXc3giFImV3Bb3d6FGo3IYiP3KXztwCLcBGAsYHQ/w466-h281/kll-comparison-x.png" width="466" /></a></div><div style="text-align: left;">In these figures, you can see that all three digests have zero error for small number of samples such that $qn \approx 10$. This is because both digests are remembering the exact sample values for the 10 smallest samples by collapsing more medial samples together except in the notable case of the smaller $t$-digest with $q=10^{-4}$ and $n=10^6$. When $qn\approx 100$, this is no longer possible and the ability to interpolate between singletons and larger centroids becomes crucial. This transition point is the worst case for the $t$-digest and for $qn > 100$, the ability to interpolate with larger centroids allows even the small $t$-digest to achieve substantially smaller errors than the KLL digest.</div><h3 style="text-align: left;">Two new papers</h3><div style="text-align: left;">The <a href="https://www.sciencedirect.com/science/article/pii/S2665963820300403" target="_blank">$t$-digest was described</a> in a recent open access paper in the journal Software Impacts. It doesn't introduce any new material, but is an open-access peer-reviewed citation if you need one.</div><div style="text-align: left;"><br /></div><div style="text-align: left;">Another article is a bit more confrontational and is as yet unpublished as far as I know but <a href="https://arxiv.org/pdf/2102.09299.pdf" rel="nofollow" target="_blank">is available in pre-print form</a>. This article concisely describes both the $t$-digest and the KLL sketch, but is notable for providing an attack on the $t$-digest accuracy in the form of a pathological distribution (appropriately named $\mathcal D_{hard}$) that is designed to be roughly as skewed as any distribution can be when represented with floating point numbers. In fact, the central half of the attack distribution is in the range $[ -10^{-154}, 10^{-154}]$ while the maximum value is $\approx 10^{308}$. For reference, the ratio of the age of the observable universe to the time light takes to transit a proton's diameter is only about 33 orders of magnitude. </div><div style="text-align: left;"><br /></div><div style="text-align: left;">Note that such an adversarially designed distribution actually <i>must</i> exist because the fixed memory guarantees of the $t$-digest are impossible to meet if you maintain relative accuracy for all distributions. Thus, the question is not <i>whether</i> there are distributions that degrade the accuracy of the $t$-digest. Rather, it is a question of whether it is plausible that any real data will ever be like these pathological cases and how gracefully $t$-digest behaves under duress. The evidence so far indicates that no real data is anything like $\mathcal D_{hard}$ and that $t$-digest is relatively well behaved even when it cannot succeed.</div><div style="text-align: left;"><br /></div><div style="text-align: left;">The $\mathcal D_{hard}$ distribution is still useful in less extreme forms for hardening the $t$-digest even if it doesn't apply much to real data. As an example, here is a graph that shows the accuracy of estimating the CDF of $\mathcal D_{hard}$ with a scale parameter of $e_{Max}=50$ which corresponds to a dynamic range of $10^{75}$. Note that the axes are on a log-odds scale which can de-emphasize errors.</div><div style="text-align: left;"><br /></div><div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-apXQw5CzOBQ/YL6mYyV4aTI/AAAAAAAAylA/JIq7Lw1OqXA1469Ju6xqZbHFsMbpDzz_ACLcBGAsYHQ/s900/dh-eMax%253D50%252Cdelta%253D200.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="800" data-original-width="900" height="355" src="https://1.bp.blogspot.com/-apXQw5CzOBQ/YL6mYyV4aTI/AAAAAAAAylA/JIq7Lw1OqXA1469Ju6xqZbHFsMbpDzz_ACLcBGAsYHQ/w400-h355/dh-eMax%253D50%252Cdelta%253D200.png" width="400" /></a></div><div style="text-align: left;"><br /></div><div style="text-align: left;"><div class="separator" style="clear: both; text-align: center;"><div class="separator" style="clear: both; text-align: center;"><div style="text-align: left;">In contrast, when the scale is pushed to the maximum level $e_{Max}=308$, we can see that the quantile estimation of a $t$-digest with compression factor of 200 is still accurate for the tails but essentially all resolution has been lost in the range from $q \in [0.1, 0.9]$. For reference, the range for the problematic region is $6 \times 10^{221}$ smaller than the total range. The range for the central half of the distribution is smaller than the total range by an amount quite a bit larger than the maximum dynamic range for 64-bit IEEE floating point numbers.</div><div style="text-align: left;"><br /></div></div></div><div class="separator" style="clear: both; text-align: center;"><div class="separator" style="clear: both; text-align: center;"><div class="separator" style="clear: both; margin-left: 1em; margin-right: 1em;"><div class="separator" style="clear: both; text-align: center;"><div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-xZK2gmoIg0Q/YL6mord3IXI/AAAAAAAAylM/9q7bVLxV1cw0i5iAf4v2qrVXM1u7jfEVACLcBGAsYHQ/s900/dh-eMax%253D308%252Cdelta%253D200.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="800" data-original-width="900" height="355" src="https://1.bp.blogspot.com/-xZK2gmoIg0Q/YL6mord3IXI/AAAAAAAAylM/9q7bVLxV1cw0i5iAf4v2qrVXM1u7jfEVACLcBGAsYHQ/w400-h355/dh-eMax%253D308%252Cdelta%253D200.png" width="400" /></a></div></div></div></div></div></div></div><p></p>Increasing $\delta$ to a value of 500 decreases the region where $t$-digest is unable to produce useful estimates to roughly the central half of the $\mathcal D_{hard}$ distribution. So, if you actually have data with this much dynamic range, you have a work-around.<h3 style="text-align: left;">Lots of code cleanups </h3><div>This release has benefited from several years to code cleanups, many due to community contributions by an army of careful readers of the code. Most of these improvements have been made to the <span style="font-family: verdana;">MergingDigest</span> to the partial neglect of the <span style="font-family: verdana;">AVLTreeDigest</span>. </div>Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com0tag:blogger.com,1999:blog-4453805095942812863.post-50292939827535580992020-10-13T13:29:00.005-07:002020-10-13T13:36:18.896-07:00What is the right number of clusters?<div style="text-align: left;"><div style="text-align: left;"><span style="color: #282829; font-family: arial;"><span style="font-size: 15px;">This is a question that comes up *all* the time when people first encounter clustering. Sadly, the information that they get about how to choose the right number of clusters is exactly wrong in a lot of the situations. The single right answer is that you should choose the parameterization of your clustering that optimizes whatever downstream use you have for clustering.</span></span></div><div><span style="color: #282829; font-family: arial;"><span style="font-size: 15px;"><br /></span></span></div><div><span style="color: #282829; font-family: arial;"><span style="font-size: 15px;">Far too often, that downstream use is simply ignored and advice like finding the knee in the residual is given as if it applied universally.</span></span></div><div><span style="font-family: arial;"><span face="-apple-system, sytem-ui, system-ui, "Segoe UI", Roboto, Oxygen-Sans, Ubuntu, Cantarell, "Helvetica Neue", sans-serif" style="background-attachment: initial; background-clip: initial; background-image: none; background-origin: initial; background-position: initial; background-repeat: initial; background-size: initial; color: #282829; font-size: 15px;">There are some really cool algorithms that apparently let the number of clusters come from the data. For instance, check this</span><span face="-apple-system, sytem-ui, system-ui, "Segoe UI", Roboto, Oxygen-Sans, Ubuntu, Cantarell, "Helvetica Neue", sans-serif" style="background-attachment: initial; background-clip: initial; background-image: none; background-origin: initial; background-position: initial; background-repeat: initial; background-size: initial; color: #282829; font-size: 15px;"> <a href="https://www.cs.cmu.edu/~kbe/dp_tutorial.pdf">tutorial on Dirichlet Processes and infinite mixture modeling</a>. </span><span face="-apple-system, sytem-ui, system-ui, "Segoe UI", Roboto, Oxygen-Sans, Ubuntu, Cantarell, "Helvetica Neue", sans-serif" style="background-attachment: initial; background-clip: initial; background-image: none; background-origin: initial; background-position: initial; background-repeat: initial; background-size: initial; color: #282829; font-size: 15px;">Or take a look at this <a href="https://icml.cc/2012/papers/291.pdf">description of DP-means</a> which applies these ideas of Bayesian modeling to the problem of clustering.</span></span></div></div><p><span style="font-family: arial;"><span face="-apple-system, sytem-ui, system-ui, "Segoe UI", Roboto, Oxygen-Sans, Ubuntu, Cantarell, "Helvetica Neue", sans-serif" style="background-attachment: initial; background-clip: initial; background-image: none; background-origin: initial; background-position: initial; background-repeat: initial; background-size: initial; color: #282829; font-size: 15px;">These are very exciting and you can use these methods for some pretty impressive practical results. For instance, <a href="https://papers.nips.cc/paper/4362-fast-and-accurate-k-means-for-large-datasets.pdf">this really fast version of $k$-means</a> </span><span face="-apple-system, sytem-ui, system-ui, "Segoe UI", Roboto, Oxygen-Sans, Ubuntu, Cantarell, "Helvetica Neue", sans-serif" style="color: #282829; font-size: 15px;">uses something akin to DP-means as an internal step to build good sketches of an empirical distribution that can give you a high quality clustering (probably).</span></span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px 0px 1em; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;">But there is a really dark side here. The fact is that you are just moving the arbitrary decision about the number of clusters into an arbitrary decision about which prior distribution to pick.</span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; direction: ltr; margin: 0px 0px 1em; overflow-wrap: anywhere; padding: 0px;"><span style="color: #282829; font-family: arial;"><span style="font-size: 15px;">The fact is, there is no universal best number of clusters for $k$-means clustering. You have to make assumptions and those assumptions may run very, very counter to reality. Whether those assumptions concern a heuristic choice of $k$ for $k$-means, or a heuristic choice of a scaling parameter $\lambda$ for DP-means, or a base distribution $G_0$ for a Dirichlet process, or even just deciding which data to ignore, you are still making a strong decision about what is best and that decision does not come from the data. Each of these decisions, together with choices like how you define distance are just different ways of encoding our prejudices.</span></span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px 0px 1em; overflow-wrap: anywhere; padding: 0px;"><span style="background-color: initial;"><span style="font-family: arial;">And it is completely necessary to encode our prejudices.</span></span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px 0px 1em; overflow-wrap: anywhere; padding: 0px;"><span style="font-family: arial;"><span style="background: none;">The real-world term for these assumptions is called “domain knowledge”. You need to know <i>why</i> y</span><span style="background: none;">ou are clustering data to understand <i>how</i> </span><span style="background: none;">to cluster it.</span></span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px 0px 1em; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;">As a classic example, the data in the figure below are all the exact same synthetic dataset viewed with different definitions of distance. Should there be one, two, three or four clusters? Or should there be 200? Each of these is a correct answer for different applications.</span></p><div class="CssComponent-sc-1oskqb9-0 QTextImage___StyledCssComponent-sc-1yi3aau-0 ilDNcJ" style="background-color: white; color: #282829; font-size: 15px;"><div class="q-box unzoomed" style="box-sizing: border-box; cursor: -webkit-zoom-in; direction: ltr; filter: blur(0px); margin-bottom: 1em;"><img class="q-image qu-display--block qu-borderRadius--small" src="https://qph.fs.quoracdn.net/main-qimg-963bb88ea2bdc1bdbb78ac7737f1e8a9" style="animation-duration: 0.001s; animation-name: insQ_100; border-radius: 3px; border: 0px none; box-sizing: border-box; color: transparent; direction: ltr; display: block; max-width: 100%;" /></div></div><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="font-family: arial;"><span style="background: none;">(see </span><span class="q-inline" style="box-sizing: border-box; direction: ltr; display: inline;"><a class="q-box qu-cursor--pointer qu-hover--textDecoration--underline" href="https://gist.github.com/tdunning/badb88043d41d916a3148c669f2fb0cd" rel="noopener nofollow" style="-webkit-tap-highlight-color: rgba(255, 255, 255, 0.6); background: none; border-radius: inherit; box-sizing: border-box; color: #195faa; cursor: pointer; direction: ltr; text-decoration-line: none; transition-duration: 180ms; transition-property: text-decoration; transition-timing-function: ease-in-out;" target="_blank" title="gist.github.com">how different definition of distance changes our view of clustering</a>)</span><span style="background: none;"> for the R script that generated these figures)</span></span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;"><br /></span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;">So let's step through some practical thoughts on clustering:</span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;"><br /></span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;">- first, think about what you want to do. Is the clustering going to make pretty pictures? Or are you using the clustering as an input to a machine learning system?</span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;"><br /></span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;">- second, think about distance and scaling. It is very unlikely that the data as you get it will be what you want. Some variables might be nearly insignificant and should be removed or scaled down. Others might be important, but measured on a very different scale from the rest of the variables you observe. There might well be other ways to improve the sense of proximity than Euclidean distance on your original data. That takes us to the next point.</span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;"><br /></span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;">- consider dimensionality reduction before clustering. Algorithms like <a href="https://umap-learn.readthedocs.io/en/latest/">UMAP</a> can vastly simply the structure of your data (if you have a decent metric). The sophistication of your projection can have the result of allowing a simpler clustering algorithm.</span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;"><br /></span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;">- try something like $k$-means or $DP$-means to see how it goes. To evaluate this, test your data against its ultimate use. If that is making pictures, look at the pictures. If it is intended to generate features for a downstream model check how much lift your clustering gives you. </span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;"><br /></span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;">As an example, suppose that our data is heavily intertwined like these two spirals.</span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;"><br /></span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;"></span></p><div class="separator" style="clear: both; text-align: center;"><a href="https://lh3.googleusercontent.com/-4u7kaP7a-28/X4YGrHuGm5I/AAAAAAAAnc8/_J_6znnD6-wXiX0axLP0StNBnkkTqM24gCLcBGAsYHQ/image.png" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><span style="font-family: arial;"><img alt="" data-original-height="878" data-original-width="1000" height="240" src="https://lh3.googleusercontent.com/-4u7kaP7a-28/X4YGrHuGm5I/AAAAAAAAnc8/_J_6znnD6-wXiX0axLP0StNBnkkTqM24gCLcBGAsYHQ/image.png" width="273" /></span></a></div><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background-color: initial; font-family: arial;">The conventional wisdom is that this kind of data is not suitable for simple clustering methods like $k$-means. But here, I have used 40 clusters (see the x marks) and you can see that the data nearest to each cluster is from a single class. That means that the distance to each of the clusters is an interesting restatement of the original data that makes things much simpler. </span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; direction: ltr; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; color: #282829; font-family: arial; font-size: 15px;"><br /></span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;">Note that this is the <i>opposite</i> of dimensionality reduction. The original data had 2 dimensions, but the restated data as distances to clusters now has 40 dimensions. Classification in this restated space is vastly simpler than in the original space because we just have to find which of the 40 values is smallest to find which color the point should be. And, indeed, that is what we see. </span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="font-family: arial;"><br /></span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;"></span></p><div class="separator" style="clear: both; text-align: center;"><div class="separator" style="clear: both; text-align: center;"><a href="https://lh3.googleusercontent.com/-ZWUAVYv9A0E/X4YIN_GMY0I/AAAAAAAAndM/DEqTWfU5qrYoyv37oUDAeE5NwmbyjSoxwCLcBGAsYHQ/image.png" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><span style="font-family: arial;"><img alt="" data-original-height="1040" data-original-width="1290" height="240" src="https://lh3.googleusercontent.com/-ZWUAVYv9A0E/X4YIN_GMY0I/AAAAAAAAndM/DEqTWfU5qrYoyv37oUDAeE5NwmbyjSoxwCLcBGAsYHQ/image.png" width="298" /></span></a></div></div><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;"></span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; direction: ltr; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="color: #282829; font-family: arial;"><span style="font-size: 15px;">Here is an ROC curve for two linear models. The brown line shows a linear classifier that is based on the original data. You can see how it barely exceeds the performance of a random classifier (the black diagonal). On the other hand, the classifier that uses the inverse distance to each of the 40 clusters is very good (you can tell because the green line goes very close to the upper left corner).</span></span></p><div><span style="font-family: arial;"><br /></span></div><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;">In fact, this model is so simple that you could hand code it by finding the nearest cluster and using the dominant color for that cluster.</span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;"><br /></span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;">This is wonderful and cool. It should be remembered, however, that you could get similarly good results by using a fairly simple neural network with enough hidden nodes. This is because the hidden nodes in a neural network function very similarly to the way that this cluster based classifier is working. This happens because the multiplication by the weights for the first layer is a dot product and this dot product is closely related to Euclidean distance.</span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;"><br /></span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;">The point here is that 40 clusters is far more than would have been recommended by most of the simple criteria that many people recommend (such as finding the knee in the total distance). That difference is because the intended use is very different.</span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;"><br /></span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none; font-family: arial;">So to mis-quote Forrest Gump, clustering is as clustering does!</span></p><p class="q-text qu-display--block" style="background-color: white; box-sizing: border-box; color: #282829; direction: ltr; font-size: 15px; margin: 0px; overflow-wrap: anywhere; padding: 0px;"><span style="background: none;"><span style="font-family: arial;"><br /><br /></span><br /></span></p>Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com0tag:blogger.com,1999:blog-4453805095942812863.post-67916375009286541762018-03-25T14:53:00.001-07:002018-03-25T14:53:05.152-07:00Getting the right answer ... much laterA whole lot of years ago, I was a poor student at the University of Colorado in Boulder. Having nothing better to do (like school work) I was wandering the halls of the upper floors of the <a href="https://www.colorado.edu/even/about-us/map-and-directions/engineering-center" target="_blank">Engineering Center</a>, and one name tag caught my eye. It said "S. Ulam".<br /><br />Being a geek and having gone to high school in Los Alamos, I knew about Stanislaw Ulam. At least, I knew a bit about the work he had done during the Manhattan project at Los Alamos. And that he was among the greats of mathematics. But it seemed totally implausible that the person behind that door could really be Stanislaw Ulam. Ulam was a creature of myth to me. Not just from another generation, but practically from another world. There was no way that it could be him. No way.<br /><br />So I had to find out.<br /><br />I knocked.<br /><br />And was invited in.<br /><br />So I asked, "Are you Stanislaw Ulam?". He seemed amused. "Yes". "Uh..." I said.<br /><br />Just as a note to anyone meeting a god who steps down onto the earth to walk among us, I suggest you have a question figured out ahead of time because you aren't going to come up with anything very good in the moment. At least I couldn't. In a similar vein, figure out your three wishes ahead of time as well in case you meet a genie.<br /><br />After a moment of hurried thought, I came up with "Could we talk?". Seriously. That is the best I could do.<br /><br />He said that would be fine, but that he had about 20 minutes of work to do first. To keep me busy, he gave me a problem to work on during that time. Not surprisingly, that problem was a really good one for assessing my mathematical sophistication.<br /><br />That sophistication was not very high (and really still isn't decades later) partly because of youth and partly because of attitude and partly because schools can't deal with students who learn quickly. My history with math was that during an explosive few weeks in the 8th grade I had gotten the bug and worked my way through all of the high school math curriculum including calculus. It was a wonderful time and didn't feel like learning so much as just remembering how things worked. But there wasn't much anything to be done beyond that because there wasn't much anybody to talk to who knew what lay past that because we lived on a military base in Germany. I could wander through books, but I didn't really have enough understanding to go much further than a bit of differential geometry and multivariate calculus. Later, when we moved to Los Alamos, there were all kinds of people I could have learned from but I was engrossed in electronics and computing and German and Russian.<br /><br />The problem that he gave me was to prove a simple fixed point theorem. Given a set $S$ and a continuous function $f: S \rightarrow S$ such that $f$ is a contraction, show that $f$ has a fixed point in $S$. A contraction is a function that maps distinct points closer together, $x_1 \ne x_2 \implies \rho(x_1, x_2) > \rho(f(x_1), f(x_2))$ where $\rho$ is some metric on $S$.<br /><br />As with all good math questions, the first response should always be clarifications to make sure that the problem is actually well understood. Even though it has been a very long time, I remember asking if I could assume that $S$ is compact and that $f$ is continuous. As I remember, I didn't actually finish a proof during the short time I had then, but I did sketch out an approach that started by noting that $\rho(x, f(x)) > \rho(f(x), f(f(x)))$. Unfortunately, I took a wrong turn at that point and thought about limits. That apparently didn't matter and I had some great talks with Professor Ulam and his friend Jan Micielski.<br /><br />One of the cool things about mathematics is that questions like this wind up in the back of your head and can pop out at any time. That just happened as I was reading a review of the latest edition of <a href="https://www.springer.com/us/book/9783642008566" target="_blank">Proofs From the Book</a>. This looks like a great book with all kinds of interesting insights to be had, but I just looked at their wonderfully short proof of the fundamental theorem of algebra. That proof goes something like this<br /><br /><ol><li>Every real-valued function $f$ over a compact set $S$ has a minimum in $S$</li><li>Every point $x$ such that a complex polynomial is non-zero, $|p(x)|>0$, has a nearby point $x'$ such that $|p(x')| < |p(x)|$</li><li>Thus, the minimum of $|p(x)|$ exists, but all points $x$ such that $|p(x)|>0$ are not the minimum. Therefore, for some point $x_0$ we have $|p(x_0)|=0$</li></ol><div>Point 2 takes a tiny bit of algebra, but overall this proof is stunningly simple and beautiful. But in the way that these things happen, after reading this proof I remembered the contraction theorem that I hadn't quite come up with so long ago and it occurred to me that this same trick would apply there. </div><div><br /></div><div>The outline is </div><div><ol><li>For any fixed point $x^*$, we have $\rho(x^*,f(x^*)) = 0$</li><li>For any other point $x$ we have $\rho(x, f(x)) > 0$</li><li>The function $d(x) = \rho(x, f(x))$ is continuous and has a minimum in $S$</li><li>For any point $x$ such that $d(x) \ne 0$, we have $d(x) > d(f(x))$ and thus $d(x)$ is not minimized at $x$</li><li>Thus, there exists some point $x^*$ where $d(x^*)=0$, that is to say $x^*$ is a fixed point.</li></ol><div>Isn't it just the way of the world that you think of the perfect answer just a little bit too late?</div></div><br /><style type="text/css">p.p1 {margin: 0.0px 0.0px 0.0px 0.0px; font: 12.0px Monaco} </style>Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com0tag:blogger.com,1999:blog-4453805095942812863.post-54401195962858075782017-06-17T15:25:00.000-07:002018-03-25T14:21:24.714-07:00Generalized Log-Likelihood Test for Poisson DistributionIn another blog, I described how a generalized log-likelihood ratio can be used to find interesting differences in counts. In the simplest cases, we have counts for some kind of observation (A and not A) under two conditions (B and not B). The question of interest is whether the the rate of A varies with whether or not condition B applies. In classical statistics, we would be talk about having a test against a null hypothesis, but in practice we want to look at lots of A's under lots of different conditions B. Testing so many situations makes it very hard to use the classical machinery well, so we have to look at the problem a bit differently.<br /><div><br /></div><div>As I mentioned in other blogs, we can still use a classically derived test known as the generalized log-likelihood ratio as a way of simply ranking different A-B combinations against each other according to how interesting they are. Even without being able to interpret the statistical score as a statistical test, we get useful results in practice.</div><div><br /></div><div>The generalized log-likelihood ratio most commonly used in these situations is derived assuming we have two binomial observations. This test can be extended to compare two multinomial conditions for independence, but this is rarely done, if only because comparing two binomials is so darned useful.</div><div><br /></div><div>With the binomial test, we look at the number of positive observations out of some total number of observations for each condition. In some situations, it is much more natural to talk about the number of positive observations not as a fraction of all observations, but as a fraction of the duration of the condition. For instance, we might talk about the number of times we noticed a particular kind of network error under different conditions. In such a case, we probably can say how long we looked for the errors under each condition, but it can be very hard to say how many observations there were without an error.</div><div><br /></div>For a Poisson distribution under two conditions $A$ and $\neg A$, we can observe a count for each of the conditions as well as the total time over which the count is taken. We can arrange our results this way:<br /><div style="text-align: left;"><br /></div><div><table align="center" border="0" cellpadding="5" cellspacing="0" style="border-collapse: collapse; text-align: left;"><tbody><tr> <td></td> <td>Count</td> <td>$$\Delta t$$</td> </tr><tr> <td><div style="text-align: right;">$$A$$</div></td> <td style="border: 1px solid #001;"><div style="text-align: center;">$$k_1$$</div></td> <td style="border: 1px solid #001;"><div style="text-align: right;">$$\,\,\,\, t_1 \,\,\,\,$$</div></td> </tr><tr> <td><div style="text-align: right;">$$\neg A$$</div></td> <td style="border: 1px solid #001;"><div style="text-align: center;">$$k_2$$</div></td> <td style="border: 1px solid #001;"><div style="text-align: right;">$$t_2$$</div></td> </tr></tbody></table></div><div style="text-align: left;"><br /></div><div style="text-align: left;"><br /></div>This is suggestive of the way that counts from two binomial observations can be arranged to look for a difference under different conditions \cite{dunning93}.<br /><br />We can investigate whether the Poisson distribution the same under both conditions using the generalized log-likelihood ratio test. Such a test uses $\lambda$, the generalized log-likelihood ratio,<br />\[<br />\lambda = \frac{<br />\max_{\theta_1 = \theta_2 = \theta_0} p(k_1 | \theta_0)p(k_2 | \theta_0)<br />}{<br />\max_{\theta_1, \theta_2} p(k_1 | \theta_1)p(k_2 | \theta_2)<br />}<br />\]<br /><div>According to Wilks\cite{wilks1938} and also later Chernoff\cite{Chernoff1954}, the quantity $-2 \log \lambda$ is asymptotically $\chi^2$ distributed with one degree of freedom.<br /><br />For the Poisson distribution,<br />\[<br />p(k | \theta, t) = \frac{(\theta t)^k e^{-\theta t}}{k!} \\<br />\log p(k|\theta, t) = k \log \theta t - \theta t - \log k!<br />\]<br />The maximum likelihood estimator $\hat\theta$ can be computed by maximizing the log probability</div><div>\[<br />\max_\theta \frac{(\theta t)^k e^{-\theta t}}{k!} = \max_\theta \log \frac{(\theta t)^k e^{-\theta t}}{k!} \\<br />\log \frac{(\theta t)^k e^{-\theta t}}{k!} = k \log \theta t - \theta t - \log k! \\ <br />\frac{\partial \log L(k | \theta, t)}{\partial \theta} = \frac{k}{ \theta} - t <br />=0 \\ <br />\hat \theta = k <br />\]<br />Returning to the log-likelihood ratio test, after some cancellation we get</div><div>\[<br />-\log \lambda =<br />k_1 \log k_1 + <br />k_2 \log k_2 <br />- k_1 \log \frac{k_1+k_2}{t_1+t_2} t_1 - k_2 \log \frac{k_1+k_2}{t_1+t_2} t_2 <br />\]</div><div>Some small rearrangement gives the following preferred form that is very reminiscent of the form most commonly to compute the log-likelihood ratio test for binomials and multinomials<br />\[<br />-2 \log \lambda = 2 \left( k_1 \log \frac{k_1}{t_1} +<br />k_2 \log \frac{k_2}{t_2} - (k_1+k_2) \log \frac{k_1+k_2}{t_1+t_2} <br />\right)<br />\]</div><div><br /><br /><br /><br /><br /><br /></div><style type="text/css">p.p1 {margin: 0.0px 0.0px 0.0px 0.0px; font: 12.0px Monaco} </style>Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com0tag:blogger.com,1999:blog-4453805095942812863.post-39213479839800721852013-09-18T15:18:00.000-07:002013-09-18T15:18:50.359-07:00How Accurate is Mahout for Summing Numbers?A question was recently posted on the Mahout mailing list suggesting that the Mahout math library was "unwashed" because it didn't use Kahan summation. My feeling is that this complaint is not founded and Mahout is considerably more washed than the original poster suggests. Here is why I think this.<br /><br />As a background, if you add up lots of numbers using a straightforward loop, you can lose precision. In the worse case the loss is \(O(n \epsilon)\), but in virtually all real examples the lossage is \(O(\epsilon \sqrt n)\). If we are summing a billion numbers, the square root is \(\approx 10^5\) so we can potentially lose 5 sig figs (out of 17 available with double precision). <br /><br />Kahan summation increases the number of floating point operations by \(4 \times\), but using a clever trick and manages to retain most of the bits that would otherwise be lost. Shewchuk summation uses divide and conquer to limit the lossage with \(O(\log n)\) storage and no increase in the number of flops.<br /><br />There are several cases to consider:<br /><br />1) online algorithms such as OnlineSummarizer. <br /><br />2) dot product and friends.<br /><br />3) general matrix decompositions<br /><br />In the first case, we can often have millions or even billions of numbers to analyze. that said, however, the input data is typically quite noisy and signal to noise ratios \(> 100\) are actually kind of rare in Mahout applications. Modified Shewchuk estimation (see below for details) could decrease summation error from a few parts in \(10^{12}\) to less than 1 part in \(10^{12}\) at minimal cost. These errors are \(10^{10}\) smaller than the noise in our data so this seems not useful.<br /><br />In the second case, we almost always are summing products of sparse vectors. Having thousands of non-zero elements is common but millions of non-zeros are quite rare. Billions of non-zeros are unheard of. This means that the errors are going to be trivial.<br /><br />In the third case, we often have dense matrices, but the sizes are typically on the order of \(100 \times 100\) or less. This makes the errors even smaller than our common dot products.<br /><br />To me, this seems to say that this isn't worth doing. I am happy to be corrected if you have counter evidence.<br /><br />Note that BLAS does naive summation and none of the Mahout operations are implemented using anything except double precision floating point.<br /><br />Here is an experiment that tests to see how big the problem really is:<br /><div><span style="font-family: Courier New, Courier, monospace;"><verbatim></verbatim></span><br /><span style="font-family: Courier New, Courier, monospace;">@Test</span><br /><span style="font-family: Courier New, Courier, monospace;">public void runKahanSum() {</span><br /><span style="font-family: Courier New, Courier, monospace;"> Random gen = RandomUtils.getRandom();</span><br /><span style="font-family: Courier New, Courier, monospace;"><br /></span><span style="font-family: Courier New, Courier, monospace;"> double ksum = 0; // Kahan sum</span><br /><span style="font-family: Courier New, Courier, monospace;"> double c = 0; // low order bits for Kahan</span><br /><span style="font-family: Courier New, Courier, monospace;"> float sum = 0; // naive sum</span><br /><span style="font-family: Courier New, Courier, monospace;"> float[] vsum = new float[16]; // 8 way decomposed sum</span><br /><span style="font-family: Courier New, Courier, monospace;"> for (int i = 0; i < 1e9; i++) {</span><br /><span style="font-family: Courier New, Courier, monospace;"> float x = (float) (2 * gen.nextDouble() - 1);</span><br /><span style="font-family: Courier New, Courier, monospace;"> double y = x - c;</span><br /><span style="font-family: Courier New, Courier, monospace;"> double t = ksum + y;</span><br /><span style="font-family: Courier New, Courier, monospace;"> c = (t - ksum) - y;</span><br /><span style="font-family: Courier New, Courier, monospace;"> ksum = t;</span><br /><span style="font-family: Courier New, Courier, monospace;"> sum += x;</span><br /><span style="font-family: Courier New, Courier, monospace;"> vsum[i % 16] += x;</span><br /><span style="font-family: Courier New, Courier, monospace;"> }</span><br /><span style="font-family: Courier New, Courier, monospace;"> // now add up the decomposed pieces</span><br /><span style="font-family: Courier New, Courier, monospace;"> double zsum = 0;</span><br /><span style="font-family: Courier New, Courier, monospace;"> for (int i = 0; i < vsum.length; i++) {</span><br /><span style="font-family: Courier New, Courier, monospace;"> zsum += vsum[i];</span><br /><span style="font-family: Courier New, Courier, monospace;"> }</span><br /><span style="font-family: Courier New, Courier, monospace;"> System.out.printf("%.4f %.4f %.4f\n", </span><br /><span style="font-family: Courier New, Courier, monospace;"> ksum, 1e6 * (sum - ksum) / ksum, </span><br /><span style="font-family: Courier New, Courier, monospace;"> 1e6 * (zsum - ksum) / ksum);</span><br /><span style="font-family: 'Courier New', Courier, monospace;">}</span><br /><span style="font-family: Courier New, Courier, monospace;"></span></div><div><div style="color: #222222; font-family: arial; font-size: small;"><br /></div><div style="color: #222222; font-family: arial; font-size: small;">A typical result here is that naive summation gives results that are accurate to within 1 part in <span style="color: black; font-family: Times; font-size: small;">\(10^{12}\)</span> 8 way summation manages \(< 0.05\) parts in \(10^{12}\) and 16 way summation is only slightly better than 8 way summation.</div><div style="color: #222222; font-family: arial; font-size: small;"><br /></div><div style="color: #222222; font-family: arial; font-size: small;">If the random numbers being summed are changed to have a mean of zero, then the relative error increases to 1.7 parts in \(10^{12}\) and 0.3 parts in \(10^{12}\), but the absolute error is much smaller.</div><div style="color: #222222; font-family: arial; font-size: small;"><br /></div><div style="color: #222222; font-family: arial; font-size: small;">Generally, it doesn't make sense to do the accumulation in float's because these operations are almost always memory channel bound rather than CPU bound. Changing to floating point arithmetic in spite of this decreases the accuracy to about 500 parts per million, 200 parts per million respectively for naive summation and 8 way summation </div><div><br /></div></div>Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com1tag:blogger.com,1999:blog-4453805095942812863.post-27074157329026797232013-04-24T20:48:00.000-07:002013-04-29T18:12:40.507-07:00Learning to Rank, in a Very Bayesian WayThe problem of ranking comments by a crowd-sourced version of "quality" is a common one on the internet.<br /><br /><a href="http://simplemlhacks.blogspot.ca/2013/04/reddits-best-comment-scoring-algorithm.html">James Neufeld suggests</a> that Bayesian Bandit algorithms can be applied to this problem. The basic idea is that you would define a stochastic quality metric whose distribution for each comment depends on the up and down votes that comment has received.<br /><br />Normal ranking algorithms try to estimate the single best value for this quality metric. Neufeld suggests that this value should be sampled from a beta distribution which models the probability that a user would mark the comment positively given that they have marked the comment at all. To present comments to a user, the metric would be sampled independently for each comment and the comments would be sorted according to the resulting scores. Different presentations would necessarily result in different orders, but as users mark comments positively or negatively, the order should converge to one where the comments presented near the top of the list have the highest probability of being marked positively.<br /><br />One very nice thing about this approach is that it doesn't waste any cycles on determining the ranking of low quality comments. Once the quality of these comments has been determined to be relatively lower than the best columns, no more learning need be done with those comments. This accelerates learning of the ranking of the best options dramatically.<br /><div><br /></div><div>This idea is interesting enough that I built a <a href="https://github.com/tdunning/bandit-ranking">quick implementation</a> which you can find on github. The main sample code there invents several hundred "comments" each with a uniformly sampled probability of getting a positive rating. The ideal behavior for ordering the comments would be to put the comment with the highest probability of getting a positive rating first and the one with the lowest probability last. The way that the program proceeds is that it picks a pageful of twenty comments to show and then proceeds to generate ratings for each of the comments on that page according to the underlying probability associated with the items displayed. The process of generating pages of comments to show and applying feedback is repeated and performance is measured.</div><div><br /></div><div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-rla34yp5NSg/UXicg71u1dI/AAAAAAAAAMA/g19wfd62Ftw/s1600/precision-ranking.png" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" height="236" src="http://2.bp.blogspot.com/-rla34yp5NSg/UXicg71u1dI/AAAAAAAAAMA/g19wfd62Ftw/s320/precision-ranking.png" width="320" /></a></div><div>Here are some results of running the program. Here we have 200 total comments, of which 20 are shown on the page that defines which comments are rated. Precision is measured here to determine how many of the best 10 comments are shown on the page. As can be seen, the system shows immediate improvement as ratings are collected. The performance rises from the initially random 10% precision and passes 50% after 30 pages of ratings.</div><div><br />As James demonstrated in his article and as others have demonstrated elsewhere, this class of algorithm is very effective for this sort of bandit problem. What is much less well known is how easily you can build a system like this.</div><h3>Try it yourself</h3><div>To run this code, you will need git, maven and java 1.7. To download the source code and compile the system, do this<br /><br /></div><div><span style="font-family: Courier New, Courier, monospace;"> $ <b>git clone git://github.com/tdunning/bandit-ranking.git</b></span></div><div><span style="font-family: Courier New, Courier, monospace;"> $ <b>cd bandit-ranking</b></span></div><div><span style="font-family: Courier New, Courier, monospace;"> $ <b>mvn package</b></span></div><div><span style="font-family: Times, 'Times New Roman', serif;"><br /></span><br /><div>This will download all dependencies of the code, compile the code and run some tests. To run the test program, do this<br /><br /><span style="font-family: Courier New, Courier, monospace;"> $ <b>java -jar target/bandit-ranking-*-with-dependencies.jar</b></span></div></div><br />The output is a thousand lines of numbers that you can drop into R, <a href="http://www.omnigroup.com/products/omnigraphsketcher/">OmniGraphSketcher</a> or even Excel to produce a plot like the one above.<br /><h3>Quick code dissection</h3><div>In <span style="font-family: Courier New, Courier, monospace;"><a href="http://tdunning.github.io/bandit-ranking/com/mapr/bandit/BanditRanking.html" target="_blank">com.mapr.bandit.BanditRanking</a>,</span> the main program for this demo, a <span style="font-family: Courier New, Courier, monospace;"><a href="http://tdunning.github.io/bandit-ranking/com/mapr/stats/bandit/BetaBayesFactory.html" target="_blank">BetaBayesFactory</a></span> is used to construct several <span style="font-family: Courier New, Courier, monospace;"><a href="http://tdunning.github.io/bandit-ranking/com/mapr/stats/bandit/BayesianBandit.html" target="_blank">BayesianBandit</a></span> objects (for average results later). This pattern can be used with other kinds of bandit factories. </div><div><br /></div><div>The BayesianBandit objects allow you to do a variety of things include sampling (<span style="font-family: Courier New, Courier, monospace;"><a href="http://tdunning.github.io/bandit-ranking/com/mapr/stats/bandit/BayesianBandit.html#sample()" target="_blank">BayesianBandit.sample</a></span>) for the current best alternative, ranking (<span style="font-family: Courier New, Courier, monospace;"><a href="http://tdunning.github.io/bandit-ranking/com/mapr/stats/bandit/BayesianBandit.html#rank(int)" target="_blank">BayesianBandit.rank</a></span>) a number of alternatives and providing training data (<span style="font-family: Courier New, Courier, monospace;"><a href="http://tdunning.github.io/bandit-ranking/com/mapr/stats/bandit/BayesianBandit.html#train(int, double)" target="_blank">BayesianBandit.train</a></span>). Sampling is used in a traditional multi-armed bandit setting such as with A/B testing. Ranking is used as it is here for getting a list of best alternatives and training is used ubiquitously for feeding back training data to the bandit.</div><div><br /></div><div>Evaluation can be done by computing precision as is done here (how many good items are in the top 20?) or by computing regret. Regret is defined as the difference between the mean payoff of the best possible choice and the mean payoff of the choice made by the bandit. For the ranking problem here, I assume that payoff of a page is the sum of the probabilities of positively rating each item on a page.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-jo8voding4c/UXlNIbDIQtI/AAAAAAAAAMk/CO70gJZcy3U/s1600/alternative-distributions.png" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" height="263" src="http://1.bp.blogspot.com/-jo8voding4c/UXlNIbDIQtI/AAAAAAAAAMk/CO70gJZcy3U/s320/alternative-distributions.png" width="320" /></a></div>The BetaBayesFactory internally uses a <a href="http://en.wikipedia.org/wiki/Beta-binomial_distribution">beta-binomial</a> distribution to model the likelihood of a positive rating for each rank. A more general alternative would be to use a <a href="http://en.wikipedia.org/wiki/Normal-gamma_distribution">gamma-normal</a> distribution. This can be done by using the GammaNormalBayesFactory instead. This extra generality comes at a cost, however, as the graph to the left shows. Here, the beta-binomial distribution results in considerably faster convergence to perfect precision than the gamma-normal. This is to be expected since the beta-binomial starts off with the assumption that we are modeling a binary random variable that can only take on values of 0 and 1. The gamma-normal distribution has to learn about this constraint itself. That extra learning costs about 50 pages of ratings. Put another way, the cumulative regret is nearly doubled by the choice of the gamma-normal distribution.<br /><br /><div class="separator" style="clear: both; text-align: center;"></div><div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-TMwZprRhUEw/UXlEAJ53NlI/AAAAAAAAAMY/J4EXBAeZsCE/s1600/no-oversampling.png" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" height="263" src="http://2.bp.blogspot.com/-TMwZprRhUEw/UXlEAJ53NlI/AAAAAAAAAMY/J4EXBAeZsCE/s320/no-oversampling.png" width="320" /></a></div><div style="text-align: start;">In order to understand what the algorithm is really doing at a high level, the graph on the right is helpful. What it shows is the number of times comments that are at different ranks are shown. What is striking here is that comments that are below the fourth page get very few trials and even on the second page, the number of impression falls precipitously relative to the first page of comments. This is what you would expect because in this experiment, it takes only a few ratings on the worst comments to know that they stand essentially no chance of being one of the best. It is this pattern of not sampling comments that don't need precise ranking that makes Bayesian Bandits so powerful.</div><div style="text-align: start;"><br /></div></div>Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com6tag:blogger.com,1999:blog-4453805095942812863.post-70211945392481948292012-10-26T12:26:00.001-07:002013-04-25T12:04:43.082-07:00References for On-line Algorithms[See also <a href="http://tdunning.blogspot.com/2013/04/learning-to-rank-in-very-bayesian-way.html">my more recent blog</a> on this talking about using bandits for ranking rated items like comments]<br /><br />I have been speaking lately about how various on-line algorithms have substantial potential for various real-time learning applications. The most notable of these algorithms are Thompson sampling for real-time handling of multi-armed bandit and contextual bandit problems and an algorithm due to Shindler, Myerson and Wang for fast $k$-means clustering of data. Since I have had a number of requests for references back to the original sources for these works, I figured a few blog posts would be a good thing. This post will describe the multi-armed bandit work and the next will describe the clustering work.<br /><h4>The Basic Problem - Multi-armed Bandits</h4>For the bandit problems, there are two basic problems to be dealt with. The first and most basic problem is that of the multi-armed bandit. In this problem, you can sample from any of a finite number of distributions and your goal is to maximize the average value of the values that you get. This can be cast into a number of practical settings in which you select which slot machine to put a quarter into, or you select which on-line ad to present to a user or you select which landing page to deliver to a user's browser should see when they visit a particular URL. It is common to simplify this case further by assuming a stationary distribution. Obviously, at least one of the distributions you are picking from has a mean equal to the large mean of any alternative. Any time you take a sample from a distribution that has a smaller mean, you fall behind the theoretical best, on average, that you could have achieved by picking from (one of) the best distributions. The degree by which you fall behind is known as the regret that you incur.<br /><br />The key to the multi-armed bandit problem is that you cannot know which distribution might have the largest mean. This means that you have to sample all of the distributions in order to estimate their means, but this implies that you have to sample from the lesser distributions in order to determine that their are, in fact, inferior. <br /><br />There are well known bounds on how well you can actually solve this problem. There are also a number of algorithms that have regret on par with these bounds or come reasonably close to these bounds. Mostly, however, these known solutions have limitations either on the number of distributions they can consider or on the complexity of the solution.<br /><br /><a href="http://www.cs.mcgill.ca/~vkules/bandits.pdf">Kuleshov and Precup</a> provide some good examples of how to compare different bandit algorithms in their paper. This <a href="https://sites.google.com/site/banditstutorial/">tutorial on bandits</a> provides a wider view of different forms of multi-armed bandit problems with a number of references.<br /><div><br /></div><div>Conspicuously missing from most lists of references, however, is all the recent work using Thompson sampling. These algorithms, which I have referred to as Bayesian Bandits, have particularly nice properties of simplicity and optimality. <a href="http://books.nips.cc/papers/files/nips24/NIPS2011_1232.pdf">Chapelle and Li</a> provide an empirical look at performance with these algorithms compared to upper confidence bound (UCB) algorithms. The last paragraph of that paper laments the lack of a theoretical analysis of these algorithms, but that lack was cured shortly in this paper by <a href="http://arxiv.org/abs/1111.1797">Agrawal and Goyal</a>. <a href="http://www.economics.uci.edu/~ivan/asmb.874.pdf">Scott</a> provided a more comprehensive view of these algorithms under the name of randomized probability matching. </div><div><br /></div><div>The idea behind Bayesian Bandits is quite simple. For each bandit, we maintain use the observations so far to build a posterior distribution for the mean of the associated payoff distribution. For binary payoffs, it is common to use a $\beta$-binomial distribution and for other cases a $\gamma$-normal distribution works well. To pick a bandit, we sample a mean for each bandit from these posterior distributions and then pick the bandit with the largest sampled mean. The new sample from that bandit gives us more data which refines the posterior distribution for that bandit. We can repeat this process as long as desired.</div><h4>Extensions to Contextual Bandits</h4><div>One of the most important characteristics of the Thompson sampling approaches (aka randomized probability matching aka Bayesian Bandits) is that they can be extended to more complex situations. One setting that I have found particularly useful involves optimizing return not just from a few bandits, but from a parameterized set of bandits that could conceivably even be infinite. The transformation from the parameters to the bandit distribution is unknown, but if we could know that, we would be able to search the parameter space to find the bandit with the highest mean payoff. </div><div><br /></div><div>This formulation is a generalization of the previous case because we can take the parameter to be an integer from $1 \ldots k$ where there are $k$ bandits and the transformation consists of the mean payoffs for each of the $k$ bandits.</div><div><br /></div><div>The algorithm in the contextual case simply consists of sampling the transformation from some posterior distribution and then solving for the parameters of the bandit that we would like to use. Some of the parameters might be fixed by the context we are working in which is where the name contextual bandits comes in.</div><div><br /></div><div>The paper by Scott alludes to this formulation, but the most approachable work on this that I know of is the paper by <a href="http://research.microsoft.com/apps/pubs/default.aspx?id=122779">Graepel, Candela, Borchert, and Herbrich</a>. In this paper, they describe the operation of AdPredictor, a system used by the Bing search engine to target ads using context.</div>Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com0tag:blogger.com,1999:blog-4453805095942812863.post-68445954076687936452012-02-08T13:57:00.001-08:002012-03-09T18:29:24.507-08:00Bayesian BanditsThe basic idea with Bayesian Bandits is to solve the problem of the explore/exploit trade-off in a multi-armed bandit by keeping a distributional estimates for the probability of payoff for each bandit, but avoiding the cost of manipulating distributions by sampling from those distributions and then proceeding on each iteration as if that were a point estimate. <br /><br />The advantage that this is that bandwidth assignments will be made according to the best estimate of the probability that each bandit could possibly be the best. This automatically causes the system to do exploration as long as it is plausible and causes the system to smoothly transition to exploitation when it becomes clear which bandit is the best. Essentially what this gives us is a Bayesian implementation of active learning.<br /><br />To illustrate this, here is a graph that shows the posterior distribution of conversion probability for two bandits where we have lots of history for one and only a little history for the other.<br /><br /><a href="http://1.bp.blogspot.com/-u4wShj4fry4/TzLnd9pDSsI/AAAAAAAAAJM/Qcg6v84EqQU/s1600/bandit-dist.png" imageanchor="1" style="clear: left; display: inline !important; float: left; margin-bottom: 1em; margin-right: 1em; text-align: center;"><img border="0" height="251" src="http://1.bp.blogspot.com/-u4wShj4fry4/TzLnd9pDSsI/AAAAAAAAAJM/Qcg6v84EqQU/s320/bandit-dist.png" width="320" /></a>You can see that the red bandit with 100 conversions out of 1000 impressions mostly like has a probability of conversion of 0.1, more or less a bit. The blue bandit with no conversions out of 10 impressions is very likely worse than the red bandit, but there is a small possibility that it is better. If we just picked the average or mode of these distributions, we would conclude that the blue bandit is worse and wouldn't give it any bandwidth without substantial mechanisms to over-ride this decision. <br /><br />On the other hand, if we estimate a conversion probability by sampling and use that sampled estimate for targeting, then we will give the blue bandit a little bandwidth and thus a chance to redeem itself.<br /><br />There are several aspects of the Bayesian Bandit algorithm that are exciting.<br /><br /><ul><li>exploration and exploitation are handled uniformly in the same framework</li><li>our internal representation encodes all of the information that we have so we don't confuse evidence of failure with lack of evidence for success</li><li>the algorithm only requires a few lines of code</li><li>the updates for the algorithm can be expressed in terms of back-propagation or stochastic gradient descent</li><li>the performance is really good.</li></ul><div><div class="separator" style="clear: both; text-align: center;"></div><div class="separator" style="clear: both; text-align: center;"></div><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-o6XtWoBwgec/T1q72aVbOiI/AAAAAAAAAKI/i0LMvw5t8y0/s1600/logN-convergence.png" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" src="http://1.bp.blogspot.com/-o6XtWoBwgec/T1q72aVbOiI/AAAAAAAAAKI/i0LMvw5t8y0/s1600/logN-convergence.png" /></a></div><span style="color: black;">As a bit of a teaser, here are a few graphs that describe how the Bayesian Bandit converges in simulated runs. The first graph shows how the average total regret for the Bayesian Bandit algorithm approaches the ideal as the number of trials increases. This experiment uses normally distributed rewards with $\sigma = 0.1$. Any algorithm that meets the optimal $O(\log n)$ convergence lower bound is said to "solve" the bandit problem.</span></div><div><div class="separator" style="clear: both; text-align: center;"></div><div class="separator" style="clear: both; text-align: center;"></div><span style="color: black;"><br /></span></div><div><span style="color: black;">The convergence here is very close to the ideal convergence rate. Note that this graph includes the convergence to optimal payoff, not the convergence knowing which is the better bandit. This is actually an interesting aspect of the problem since the algorithm will converge almost instantly for cases where the conversion probabilities are highly disparate which will make the payoff converge quickly. For cases where the conversion probabilities are nearly the same, it will take a long time for the algorithm to determine which is the better bandit, but exploration is not expensive in such a case so the convergence to near-optimal payoff will be even faster than the case where the conversion rates are very different.</span></div><div><br /></div><div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/-FZ8WR-4_WMI/TzLulisbi_I/AAAAAAAAAJc/KddmXvs0M9k/s1600/convergence.png" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" height="251" src="http://3.bp.blogspot.com/-FZ8WR-4_WMI/TzLulisbi_I/AAAAAAAAAJc/KddmXvs0M9k/s320/convergence.png" width="320" /></a></div><div><div><span style="color: black;">For example, here is a graph of the probability of picking the better bandit where the conversion rates are nearly the same. As you can see, it takes quite a while for the algorithm to split these two options. The average payoff, however, only changes from 0.11 to 0.12 during this entire convergence and it has already reached 0.118 by the time it is 20% into the process so the cost of a long experiment is not that high.</span></div></div><div><span style="color: black;"><br /></span></div><div><span style="color: black;">Sample code for the Bayesian bandit is available at </span><a href="https://github.com/tdunning/storm-counts">https://github.com/tdunning/storm-counts</a>.</div>Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com1tag:blogger.com,1999:blog-4453805095942812863.post-73359204887254202532011-06-22T14:06:00.000-07:002011-06-22T14:06:28.390-07:00Buzzwords Keynote - Conclusion<div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">In the end, it is up to us to make things better. We need a way for non-Apache entities to interact with Apache productively. If we can't do that, then it is quite possible that all of the momentum and excitement that Hadoop now has will be lost. </div><div><br /></div><br /><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: 'Helvetica Neue', Arial, Helvetica, sans-serif; font-size: large;">Conclusion</span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">The key is that we now have an eco-system, not just a community. We can make it work. Or we can elect to let it not work. Not working is the default state. We have to take positive action to avoid the default.</div></div><div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">Apache can stay a strong voice for business friendly open source by remaining Apache. Trying to make Apache broad enough to include all of the players in Hadoop and Hadoop derivatives will simply debase the voice of Apache into the average of opposing viewpoints, i.e. into nothing.</div></div><div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div></div><div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">There <i>are,</i> however, many other players who are not part of Apache and who should probably not be part of Apache. There needs to be a way for these others to engage with the Apache viewpoint. It can't just be on the level of individuals from Apache trying to informally spread the Apache way even though that is critical to have. It is likely to require a venue in which corporate entities can deal with something comparable to themselves. A good analogy is how Mozilla participation in W3C has made the web a better place.</div></div><div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">But we can make our eco-system work. It isn’t what it was and it never will be again.</div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">But it can be astonishing</div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">Let's make it so.</div></div><div><br /></div>Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com0tag:blogger.com,1999:blog-4453805095942812863.post-47223295461962843992011-06-22T14:04:00.000-07:002011-06-22T14:04:46.587-07:00Buzzwords Keynote - Part 3<div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: Times, 'Times New Roman', serif;">In the third part of my talk, I talked a bit about where Hadoop has come from and where it is going. Importantly, this involves a choice about where Hadoop and the related companies products and individuals might be able to take things.</span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: Times, 'Times New Roman', serif;"><br /></span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: Times, 'Times New Roman', serif;"><br class="Apple-interchange-newline" /></span><span class="Apple-style-span" style="font-family: 'Helvetica Neue', Arial, Helvetica, sans-serif; font-size: large;">Where we are and how we got here</span></div><div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">My second section described the rough state of the Hadoop eco-system is a slightly provocative way. In particular, I described a time when I was on a British train and in partial compensation for delays the operators announced that "free beer would be on sale in the galley car". Free beer for sale is a wonderful analogy for the recent state of Hadoop and related software.</div></div><div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">That said, there are serious problems brewing. The current world of Hadoop is largely based on the assumption that the current community is all that there is. This is a problem, however, because the current (Apache-based) community presumes interaction by individuals with a relatively common agenda. More and more, however, the presence of a fundable business opportunity means that this happy world of individuals building software for the greater good has been invaded by non-human, non-individual corporations. Corporations can't share the same agenda as the individuals involved in Apache and Apache is constitutively unable to allow corporate entities as members.</div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">This means that the current community can no longer be the current world. What we now have is not just a community with shared values but is now an eco-system with different kinds of entities, multiple agendas, direct competition and conflicting goals. The Apache community is one piece of this eco-system.</div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: 'Helvetica Neue', Arial, Helvetica, sans-serif; font-size: large;">Our choice of roads</span></div></div><div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">Much as Dante once described his own situation, Hadoop now finds itself in the middle of the road of its life in a dark wood. The members of the Apache community have a large voice in the future of Hadoop and related software.</div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">As a darker option, the community can pretend that the eco-system that now exists of human and corporate participants is really a community. If so, it is likely that the recent problems in moving Hadoop forward will continue and even get worse. Commit wars and factionalization are likely to increase as corporate entities, denied a direct voice in Apache affairs, will tend to gain influence indirectly. Paralysis in development will stall forward progress of Hadoop itself leading to death by a thousand forks. Such a dark world would let alternative frameworks such as Azure to gain footholds and possibly to dominate.</div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">In this brighter alternative future, I think that there are ways to create a larger forum in which corporate voices can be heard in their true form rather than via conflicts of interest. In this scenario, Apache would be stronger because it really can be a strong voice of the open source community. Rather than being the average of conflicting views, Apache would be free to express the shared values of open source developers. Corporations would be able to express their goals, some shared, some not in a more direct form and would not need so much to pull the strings of Apache committers. Importantly, I would hope that Hadoop could become something analogous to a reference implementation and that commercial products derived from Hadoop would have a good way to honor their lineage without finding it difficult to differentiate themselves from the original. Hopefully in this world innovation would be welcomed, but users would be able to get a more predictable experience because they would be able to pick products offering whatever innovation rate/stability trade-off that they desire. Importantly, there would be many winners in such a world since different players would measure success in different terms.</div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">We have a key task ahead of us to define just what kind of eco-system we want. It can be mercenary and driven entirely be corporate goals. This could easily happen if Apache doesn't somehow facilitate the creation of a forum for eco-system discussion. In such an eco-system, it is to be expected that the companies that have shown a strong talent at dominating standards processes and competing in often unethical ways will dominate. My first thought when I imagine such a company is Microsoft, but that is largely based on having been on the receiving end of their business practices. I have no illusions that talent for that kind of work is exclusively found in Redmond.</div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">In my talk, I proposed some colorful cosmological metaphors for possible worlds, but the key question is how we can build a way for different kinds of entities to talk. It is important to recognize different values and viewpoints. Apache members need to understand that not everything is based on individual action, nor do corporation hold the same values. Companies need to take a strong stance to recognize the incredible debt owed to the Apache community for creating the opportunities we all see.</div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">If we can do this, then Hadoop (and off-spring) really does have a potential to dominate business computing.</div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div></div>Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com0tag:blogger.com,1999:blog-4453805095942812863.post-10970076483472876652011-06-22T14:02:00.000-07:002011-06-22T14:02:31.078-07:00Buzzwords Keynote - Part 2<div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: Times, 'Times New Roman', serif;">In the first part of the talk, I made the case that Apache Hadoop has lots of head-room in terms of performance. This translates into lots of opportunity both for open source developers to make Hadoop itself better, but also for companies to build products that derive from Hadoop but improve it in various ways.</span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: 'Helvetica Neue', Arial, Helvetica, sans-serif; font-size: large;"><br class="Apple-interchange-newline" />The $S$ score</span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">In honor of Steve Jobs whose highest praise is reputedly to say "that doesn't suck", I proposed an $S$ score whose highest score is zero, but for all real systems is always negative. For a batch, data processing system like Hadoop, I proposed that a good definition of $S$ was the log base 10 of the ratio of the actual performance to the performance implied by hardware limits.</div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">Not suprisingly, the overall score for Hadoop comes out to be somewhere between -5 to -2 depending on desired workload (i.e. Hadoop runs programs somewhere between 100 and 100,000 times slower than the hardware would allow). For some aspects, Hadoop's $S$ score can be as good as $-0.5$ but generally there are multiple choke-points and some of these are additive. This is hardly news and isn't even a mark of discredit to Hadoop since the developers of Hadoop have always prized getting things to work and to work at scale above getting things to work within an iota of the best the hardware can do at a particular scale. Another factor that drives $S$ down for Hadoop is the fact that the hardware we use has changed dramatically over the 6-7 year life of Hadoop.</div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">In defining the value of $S$ for current Hadoop versions, I don't mean to include algorithm changes. Michael Stonebraker has become a bit famous for running down Hadoop for not doing database-like things with database-like algorithms, but I would like to stick to the question of how fast Hadoop could do what Hadoop is normally and currently used to do.</div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">The key conclusion is that having such a low $S$ combined with high demand for Hadoop-like computation represents a lot of opportunity. This opportunity involves opportunities for the open source community to make things better. It also represents opportunities for commercial companies to make money. The latter kind of opportunity is what is going to shake up the currently cozy Hadoop community the most.</div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div>Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com0tag:blogger.com,1999:blog-4453805095942812863.post-17615424898158618732011-06-22T14:00:00.000-07:002011-06-22T14:09:08.069-07:00Buzzwords Keynote ... blog editionThere has been a bit of demand for an expanded version of my Buzzwords keynote from a few weeks ago. This demand has been increased by a particular unfortunate mis-quote in a tweet that suggested that I thought that there was a need for a new organization "to supersede Apache". Of course, I suggested nothing of the sort so it is a good idea to walk through the ideas that I presented. The buzz words site has the <a href="http://berlinbuzzwords.de/content/keynotes-online">video of my talk</a> and <a href="http://berlinbuzzwords.de/sites/berlinbuzzwords.de/files/buzz-words-ted-dunning.pdf">pdf of my slides</a> in case you want to follow along.<br /><br />The talk was divided into several sections. The first one proposed the uncontroversial thesis that Hadoop performs at a level far below the potential offered by modern hardware. A second section pointed out difficulties with the current social structure surrounding the development of Hadoop and related software. I then examined what I see as possible futures while describing how I think we will be choosing between these alternative futures. I will post each section in a separate blog entry.<br /><br /><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">As I spoke, I encouraged the audience to tweet using the hash-tag #bbuzz and to keep communal notes on a <a href="http://tinyurl.com/buzzwords-ted-dunning">shared Google document</a>. The tweets are a bit hard to find as befits ephemeral media, but the shared notes are still accessible.<br /><br /></div>Sections:<br /><br /><a href="http://tdunning.blogspot.com/2011/06/buzzwords-keynote-part-2.html">The S Score</a><br /><a href="http://tdunning.blogspot.com/2011/06/buzzwords-keynote-part-3.html">Possible Futures</a><br /><a href="http://tdunning.blogspot.com/2011/06/buzzwords-keynote-conclusion.html">Conclusion</a><br /><br /><br /><div><div style="direction: ltr; margin-bottom: 0pt; margin-left: .38in; margin-top: 7.68pt; mso-line-break-override: none; punctuation-wrap: hanging; text-align: left; text-indent: -.38in; unicode-bidi: embed; word-break: normal;"></div></div>Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com0tag:blogger.com,1999:blog-4453805095942812863.post-28259056207510379202011-06-20T21:06:00.000-07:002011-06-20T21:06:13.913-07:00Buzzwords WrapupWell, Buzzwords is over and my primary conclusion is that I wish I had come last year as well as this year. Isabel and Simon are really making Buzzwords into a major open source conference and with the demise of the European ApacheCon, Buzzwords is probably the the first or second most important open source conference in Europe. If you only can choose one, I would strongly recommend geeking in Berlin. It isn't just the conference; there are bunches of related events such as informal dinners, bar camps and hackathons. Since Buzzwords makes such a strong effort to include North American participants you may even have a better chance of connecting globally by going to Europe than going to a conference in the US or Canada.<br /><br />The conference itself consisted of two days of scheduled events anchored by keynotes each day. Doug Cutting gave the first keynote and covered a lot of the history and current state of Hadoop. As always, his talk was very well done and contained quite a bit of technical information which is refreshing in a keynote. I gave the second keynote and talked a bit about the state and future of Hadoop, related Apache projects and the burgeoning commercial marketplace. Some of what I said stirred up a bit of talk, which is good since my primary thesis that we aren't talking enough about how the world of Hadoop and related software is rapidly changing in ways that aren't well recognized. Stay tuned here for a blog edition of my talk.<br /><br />There were quite a few excellent technical talks as well. Among the scheduled talks, Jonathan Gray gave a talk which his usual and customary dose of excellent technical information about how Facebook is using Hbase. A notable moment came when he was asked about the state of Cassandra at Facebook. Check out the upcoming video for details on his answer.<br /><br />Dawid Weiss gave an excellent talk on finite state automata and the difference between deterministic and non-deterministic variants. The only defect I could see in his presentation was that we couldn't see the eagles on the coins. Based on the fact that the room was packed (I sat in the aisle on the floor) and the very eager audience questions, I would say that there is a surprisingly strong market place for information on foundational algorithms like finite state transformers. <br /><br />The lightning talks at the end of day two also had some gems. Thomas Hall's northern accent blended charmingly with the frank assessment of some of his experiences with certain technical approaches. I can't possibly convey the tone and content so, yet again, you will need to refer to his slides and the video on the conference web-site.<br /><br />Frank Scholten also had a lightning talk that contained a very nice walk-through of Mahout document clustering. What he showed is a work in progress, but already what he has provides a highly requested set of recipes to illustrate a lot of the software in Mahout.<br /><br />Outside of the conference there was an (excellent) barcamp run by Nick Burch. I think I learned as much about how to run a barcamp by watching him as anybody did from any of the technical discussions and the technical discussions were pretty excellent. <br /><br />I have to say that if you want to see me next year in early June, there is a high likelihood that you will have to be in Berlin to do it.<br /><br />See http://berlinbuzzwords.de/wiki/linkstoslides to get slides from talks.Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com0tag:blogger.com,1999:blog-4453805095942812863.post-31804673912105020582011-06-08T16:04:00.000-07:002011-06-08T16:04:25.799-07:00The Best Illustration of a probability DistributionDescribing a probability distribution in the abstract to a novice is often difficult.<br /><br /><a href="http://4.bp.blogspot.com/-wc8qEOoyL1s/Te__Z13wf6I/AAAAAAAAAHc/DT4JjT9IKsU/s1600/pdf-small.jpg" imageanchor="1" style="clear: right; display: inline !important; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" src="http://4.bp.blogspot.com/-wc8qEOoyL1s/Te__Z13wf6I/AAAAAAAAAHc/DT4JjT9IKsU/s1600/pdf-small.jpg" /></a>Here it is in concrete form. Note how some keys are more worn than others. This is in Germany so that the ground floor is labeled "0". Note the wear on the "door close" button!Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com4tag:blogger.com,1999:blog-4453805095942812863.post-32463476200557342682011-06-08T16:00:00.000-07:002011-06-08T23:17:47.086-07:00Visit to DIMA at Technische UniversitaetI had a great visit today at the <a href="http://www.stratosphere.eu/">DIMA laboratory at TU in Berlin</a>. They are working on an interesting system called Stratosphere which provides an interesting generalization generalization of map-reduce. Of particular interest is the run-time flexibility for adapting how the flow partitions or transfers data.<br /><br />They accomplish this by having a lower level abstraction layer that supports a larger repertoire of basic options beyond just map and reduce. These operations include match, cross product and co-group. Having a wider range of operations and retaining some additional flow information at that level allows them to do on-the-fly selection of the detailed algorithm for different operations based on the statistics of the data and the properties of the user-supplied functions.<br /><br />Here's a pic of me answering questions about startups and log-likelihood ratio tests.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://4.bp.blogspot.com/-xDvcKsIjVEE/Te_--oK2ndI/AAAAAAAAAHY/uNPMuLYOTLw/s1600/dima-visit-small.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="191" src="http://4.bp.blogspot.com/-xDvcKsIjVEE/Te_--oK2ndI/AAAAAAAAAHY/uNPMuLYOTLw/s320/dima-visit-small.jpg" width="320" /></a></div>Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com2tag:blogger.com,1999:blog-4453805095942812863.post-36195414493362872452011-06-01T18:02:00.000-07:002011-06-01T22:09:55.635-07:00A Tour of the Multi-update For ZookeeperI have recently been working on some new features for Apache Zookeeper that represent probably the largest change in the Zookeeper client API in several years. Other recent changes such as observers and read-only clusters have changed the server-side semantics, but the client API is nearly unchanged since the original public release of Zookeeper several years ago. The JIRA covering the overall issue is <a href="https://issues.apache.org/jira/browse/ZOOKEEPER-965">Zookeeper-965</a> and you can look at the code as it is being refined into committable state in my fork of Zookeeper that I keep <a href="https://github.com/tdunning/zookeeper">on github</a>.<br /><br />(<a href="http://mapr.com/?p=83&option=com_wordpress&Itemid=134">related announcement</a>)<br /><br /><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: 'Helvetica Neue', Arial, Helvetica, sans-serif; font-size: large;">The Problem</span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div>Almost from the beginning of the Zookeeper project, there have been repeated questions on the mailing about how to update several nodes at once. The answer has always been to consolidate the structures that you need to update atomically into the contents of a single znode and then update those contents atomically using standard methods. This update problem is often exacerbated by a desire to use ephemeral nodes so that state disappears correctly when a service crashes.<br /><br />This is a pretty reasonable answer and it handle many important cases fairly well. For instance, in the common case of servers that are assigned certain roles and then in turn advertise which roles they have successfully taken on, this pattern can be used by giving each server a single file of assignments and a single file of advertised capabilities. Each of these files can be ephemerally created by the server so they will both vanish when the server disappears. Assignments can be added atomically and only the server ever modifies the advertised roles. Zookeeper works quite well in these cases.<br /><br />Keeping track of a graph full of nodes with directional connections to other nodes is a good example of where Zookeeper's update model is not so nice. In this case, nodes have a list of connections to other nodes and a list of connections from other nodes. If any node has a list of outgoing connections that are not matched by the corresponding list of incoming connections on other nodes, then the graph is referentially inconsistent. Keeping the graph consistent under changes is not possible with normal Zookeeper updates unless you store the entire graph in a single file.<br /><br /><span class="Apple-style-span" style="font-family: 'Helvetica Neue', Arial, Helvetica, sans-serif; font-size: large;">Lazy Operations</span><br /><br />The new <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">multi() </span>method allows such a data structure to be updated with strong guarantees of consistency. This is done by currying the normal database mutation operators and then passing all of the resulting closures to Zookeeper at once for execution. Obviously this requires a bit of syntactic sugar in languages like Java and C which don't like to express closures directly.<br /><br />The way that this works is that there is a new class call <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Op</span>. There are sub-classes of <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Op</span> called <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Op.Create</span>, <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Op.SetData</span>, <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Op.Check</span> and <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Op.Delete</span> that correspond to the operations normally invoked by calling the similarly named methods on a ZooKeeper object. In essence, these sub-classes of Op represent reified methods or closures that can be frozen at one point in time and then executed later. These sub-classes are instantiated using factory methods on <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Op </span>class. The names and signatures of these factory methods mimic the corresponding methods on ZooKeeper.<br /><br />Once you have all of the operations you would like to perform, you can execute them all in a single transaction using <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">ZooKeeper.multi().</span> This will either execute all of the <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Op</span>'s or none of them.<br /><br /><br /><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: 'Helvetica Neue', Arial, Helvetica, sans-serif; font-size: large;">An Example</span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><a href="http://1.bp.blogspot.com/-Zfjultl0Hxs/TeRllzvRoxI/AAAAAAAAAHQ/3aYjZhX1P-8/s1600/hyper-cube.png" imageanchor="1" style="clear: right; display: inline !important; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" height="312" src="http://1.bp.blogspot.com/-Zfjultl0Hxs/TeRllzvRoxI/AAAAAAAAAHQ/3aYjZhX1P-8/s320/hyper-cube.png" width="320" /></a>I have placed an example program for doing a <a href="https://github.com/tdunning/graph-demo">graph-based computation</a> over on github. This program builds a graph consisting of 32 nodes in the form of a 5-dimensional hyper-cube and then uses a numerical technique called over-relaxation to solve for the voltages that would result if all the links in the hyper-cube were replaced by unit resistors. A picture of the graph is just to the right. In this graph, the voltage for node 0x1F is labeled as $V_5$ and the voltage for node 0x0 is labeled as $V_0$. There are lots of connections and solving this problem analytically is difficult unless you twig to the trick of using symmetry. Real-world networks generally don't exhibit such congenial properties and can be nearly impossible to solve analytically.</div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">Normally, of course, we wouldn't actually use Zookeeper for this computation, but the point here isn't so much a realistic computation as a test-bed for the new multi() method and a demonstration of how the closure based style associated with multi-ops makes code easier to read and understand.</div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: 'Helvetica Neue', Arial, Helvetica, sans-serif; font-size: large;">The Code</span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: 'Helvetica Neue', Arial, Helvetica, sans-serif; font-size: large;"></span></div><br /><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"></div><div style="font-family: Times; margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: Times; font-size: small;">To read the code, start with the class </span><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">ZGraphTest</span><span class="Apple-style-span" style="font-family: Times; font-size: small;">. This is a JUnit test that drives all the rest of the code. In this test, a graph is constructed (named </span><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace; font-size: small;">zg</span> in the code), nodes are connected in the pattern of a hyper-cube which means that nodes are labeled with integers from 0 to 31 and a node $n_1$ is connected to another node $n_2$ if $n_1$ and $n_2$ differ in exactly one bit.</div><div style="font-family: Times; margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: Times;">In each iteration of the code, the </span><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">ZGraph.reweight()</span> method of the graph is called on a randomly selected node. This sets the value at that node to be the average of the values at the neighbors of that node. This computation converges geometrically to the correct value with errors decreasing by a factor of 5-10 with every 1000 iterations. As the actual computation proceeds the error versus the theoretically known correct values is shown every 1000 iterations and you can see the errors go to zero with 6 significant figures at about 7000 iterations.</div><div style="font-family: Times; margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: Times;">Internally, a </span><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">ZGraph</span> keeps a connection to Zookeeper and the name of the root directory for all of the nodes in the graph. Methods like <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">connect()</span>, <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reweight()</span> and <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">update()</span> all work by defining lazy update or version check operations on <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Node</span> objects and then executing all of the update or version checks at one time in a transaction. For instance, in the <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">reweight()</span> method, this code gets the weight of all of the neighbors of node g1:</div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"> double mean = 0;</span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"> Set<integer> neighbors = g1.getIn();</integer></span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"> List<op> ops = Lists.newArrayList();</op></span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"> for (Integer neighbor : neighbors) {</span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"> Node n = getNode(neighbor);</span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"> ops.add(n.checkOp(root));</span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"> mean += n.getWeight();</span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"> }</span></div><div style="font-family: Times;"><br /></div><div style="font-family: Times;">As a side effect, we also collect a list of version check operations into the list ops. Then in this code we add one more operation to set the weight on g1 and add the update operation to the list ops:</div><div><div><br /></div><div><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"> // set up the update to the node itself</span></div><div><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"> g1.setWeight(mean / neighbors.size());</span></div><div><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"> ops.add(g1.updateOp(root));</span></div><div style="font-family: Times;"><br /></div><div><span class="Apple-style-span" style="font-family: Times;">Finally, </span><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">multi()</span> is called to actually do all of the version checks and updates that we have collected,</div></div><div><div><br /></div><div><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"> zk.multi(ops);</span></div><div style="font-family: Times;"><br /></div><div><span class="Apple-style-span" style="font-family: Times;">The essential point here is how the list of operations was collected a bit at a time and then executed all at once. The versions of the nodes that were inputs into the operation were collected in the form of check operations. When those check operations are executed along with the update of </span><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">g1</span>, we guarantee that <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">g1</span>'s new weight will accurately reflect the average of its neighbors and if somebody changes the value of any of the neighbors in between the time that we read the value and the time we actually do the update, we will run the update again.</div></div><div><br /></div><div>Similarly, and closer to the idea of maintaining referential integrity, the <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">ZGraph.connect()</span> collects update operations for the two nodes being connected and executes both updates together to avoid the possibility that anyone would see a situation where one node connects to a second but the second doesn't have the reverse connection. This code does the job,</div><div><div><br /></div><div><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"> Node g1 = Node.readNode(zk, root, id1);</span></div><div><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"> Node g2 = Node.readNode(zk, root, id2);</span></div><div><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"> g1.connectTo(id2);</span></div><div><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"> g2.connectFrom(id1);</span></div><div><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"> zk.multi(Arrays.asList(g1.updateOp(root), g2.updateOp(root)), results);</span></div><div><br /></div></div><div style="font-family: Times; margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">Again, the closure passing style makes this operation very easy.</div><div style="font-family: Times; margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: Times;">One additional point to be noted is that having each </span><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Node</span><span class="Apple-style-span" style="font-family: Times;"> return closures for updates that can be executed in any context the caller would like has the effect of removing all direct Zookeeper operations from the </span><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Node</span> other than for reading or creating a node. It also removes all references to serialization of a <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Node</span> from any outside code. This simplifies both the caller and the <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Node</span> because the <span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">Node</span> doesn't have to implement all plausible combinations of multiple updates and the caller doesn't have to worry about any aspect of serialization and can focus on the appropriate task of arranging the transactional semantics for updates.</div><div style="font-family: Times; margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><br /></div><div style="font-family: Times; font-size: medium; margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"></div><br /><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: 'Helvetica Neue', Arial, Helvetica, sans-serif; font-size: large;">Credit Where Credit is Due</span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-size: large;"></span></div><div style="font-family: Times; font-size: medium; margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-size: large;"><br /></span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span"><span class="Apple-style-span" style="font-family: inherit;">Of course, this new Zookeeper capability wouldn't be possible if the Apache Zookeeper project team hadn't built Zookeeper in a very flexible way in the first place. Kudos to Ben Reed and Mahadev and the other early committers on the project.</span></span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: inherit;"><br /></span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: inherit;">Also, the actual <a href="https://issues.apache.org/jira/browse/ZOOKEEPER-965">Zookeeper-965</a> project would have stalled out if Marshall McMullen and Camille Fournier hadn't stepped in with a bit of extra momentum. I had completed the wire protocols and all of the client side work, but Marshall did all of server side work and Camille provided really excellent advice about how the changes should be done.</span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-family: inherit;"><br /></span></div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">The final credit goes to MapR Technologies for supporting open source by allowing capabilities like multi to be contributed back.</div><div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;"><span class="Apple-style-span" style="font-size: large;"><br /></span></div>Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com0tag:blogger.com,1999:blog-4453805095942812863.post-51562140653094911002011-05-29T11:51:00.000-07:002011-05-31T15:30:35.804-07:00Online algorithms for boxcar-ish moving averagesOne problem with exponentially weighted moving averages is that the weight that older samples decays sharply even for very recent samples. <br /><div class="separator" style="clear: both; text-align: center;"><a href="http://4.bp.blogspot.com/-hL8j1kS5toU/TY0kbm7e8RI/AAAAAAAAAG4/R_8_VdbTaT8/s1600/graph1.png" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" height="320" src="http://4.bp.blogspot.com/-hL8j1kS5toU/TY0kbm7e8RI/AAAAAAAAAG4/R_8_VdbTaT8/s320/graph1.png" width="320" /></a></div>The impulse response of such an average shows this clearly. In the graph to the right, the red line is the impulse response of the exponential weighted average is shown by the red line. The impulse response of a different kind of moving average derived by John Canny in the early 80's is shown by the black line.<br /><br />The Canny filter puts higher weight on the events in the recent past which makes it preferable when you don't want to forget things right away, but do want to forget them before long and also want an on-line algorithm.<br /><br />The cool thing about the Canny filter is that it is only twice as much work as a simple exponential moving average. The idea is that if you take the difference between two exponentially weighted averages with different time constants, you can engineer things to give you an impulse response of the sort that you would like. <br /><br />The formula for such a difference looks like this<br />\begin{eqnarray*}<br />w(x) &=& k_1 e^{-x/\alpha} - k_2 e^{-x/\beta}<br />\end{eqnarray*}<br />Here $k_1$ and $k_2$ scale the two component filters in magnitude and $\alpha$ and $\beta$ are the time constants for the filters.<br /><br />It is nice to set the magnitude of the filter at delay $0$ to be exactly 1. We can use this to get a value for $k_2$ in terms of $k_1$<br />\begin{eqnarray*}<br />w(0) &=& k_1 - k_2 = 1 \\<br />k_2 &=& k_1 -1 <br />\end{eqnarray*}<br />Similarly, we can constrain the slope of the impulse response to be $0$ at delay $0$. This gives us $\beta$ in terms $\alpha$<br />\begin{eqnarray*}<br />w'(0) &=& {k_1 \over \alpha} - {k_2 \over \beta} = 0\\<br />{k_1 \over \alpha} &=& {k_1-1 \over \beta} \\<br />\beta &=& \alpha \frac{ k_1-1} { k_1}<br />\end{eqnarray*}<br />The final result is this impulse response<br />\begin{eqnarray*}<br />w(x) = k \exp \left(-{x \over \alpha}\right)-(k-1) \exp\left(-{k x\over \alpha (k-1)}\right)<br />\end{eqnarray*}<br />We can do a bit better if we note that as $k \to \infty$, the shape of the impulse quickly converges to a constant shape with mean of $w(x) \to \frac{3a}{2}$ and a total volume of $2a$.Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com0tag:blogger.com,1999:blog-4453805095942812863.post-73383004173696560492011-03-24T12:24:00.000-07:002011-04-08T09:11:41.411-07:00Exponentially weighted averaging for ratesIn a previous post, I talked about how to produce exponentially weighted averages of time-embedded values. This is fine for cases where you really want averages, but it doesn't work for rates. Happily, the technique I presented for simple averages can be easily extended to work for rates as well.<br /><br />To recap a bit, an on-line algorithm for computing a simple average works by accumulating the number of observations and the sum of the observations. There is no mention of time. For each observation $(x, t)$ we update our state $(s, w)$ according to<br />\begin{eqnarray*}<br />s_{n} &=& s_{n-1} + x \\<br />w_{n} &=& w_{n-1} + 1<br />\end{eqnarray*}<br />The exponentially weighted average approach that I mentioned before uses the time of each observation to discount the previous state based on the time of the previous observation. Thus, the state $(s,w,t)$ is updated according to<br />\begin{eqnarray*}<br />\pi &=& e^{-(t-t_{n-1})/\alpha} \\<br />s_{n} &=& \pi s_{n-1} + x \\<br />w_{n} &=& \pi w_{n-1} + 1 \\<br />t_{n} &=& t<br />\end{eqnarray*}<br />If we were to compute simple rates without discounting, then instead of interpreting $w$ as a count, we would interpret it as a time interval. Thus we would update state $(s,w,t)$ according to with an observation $(x, t)$ according to <br />\begin{eqnarray*}<br />s_{n} &=& s_{n-1} + x \\<br />w_{n} &=& w_{n-1} + (t - t_n) \\<br />t_n &=& t<br />\end{eqnarray*}<br />By analogy with the discounted version of the simple average, we can produce an exponentially weighted rate average by updating the state according to this<br />\begin{eqnarray*}<br />\pi &=& e^{-(t-t_{n-1})/\alpha} \\<br />s_{n} &=& \pi s_{n-1} + x \\<br />w_{n} &=& \pi w_{n-1} + (t-t_n) \\<br />t_{n} &=& t<br />\end{eqnarray*}<br />This approach has a defect in that each data pont should really be considered to be a report of events spread out in an uncertain way over the time since the last report. As such, the interval and the events should probably both be discounted as if they had been reported uniformly over the period in question. In practice, this correction does not matter much since the averaging time constant $\alpha$ is typically large compared to the average reporting interval. Another consideration comes up when multiple sources are reporting concurrently. If we want to do proper discounting of each observed number, then we have to keep track of the time since last report for each of the sources. Since this probably won't matter and since it considerably complicates the implementation, I would rather not do it.<br /><br />See https://issues.apache.org/jira/browse/MAHOUT-634 for code. This will be in Mahout before long.Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com3tag:blogger.com,1999:blog-4453805095942812863.post-36225663158597075422011-03-24T10:25:00.000-07:002011-03-24T10:25:19.507-07:00Update on exponential time-embedded averagesI will be adding this code to Mahout shortly.<br /><br />See https://issues.apache.org/jira/browse/MAHOUT-634 for status.<br /><br />Also, if you are measuring rates, then it is common for rates to be reported from multiple sources independently. Such an average can be computed pretty easily using this same framework if the sources report often relative to the averaging time constant. This simple implementation just attributes each reported count as if they occurred in the interval since the most recent report from any reporter. If the time constant is relatively long, this can work out reasonably well as long as we are careful.<br /><br />If reporting intervals are longer, then the averaging is a bit trickier because we really would like to attribute the reported counts over the entire interval from the last report from the same source. This means that we have to discount some of the counts because they are effectively kind of old.<br /><br />More details shortly.Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com0tag:blogger.com,1999:blog-4453805095942812863.post-84268918908858350652011-03-22T17:15:00.003-07:002012-03-16T17:23:22.321-07:00Exponential weighted averages with irregular samplingTed Yu on the hbase list wants to compute an exponential weighted average for the number of queries per second or other statistics that characterize the performance or state of an hbase system.<br /><br />The trick is that the measurements are only available at irregular intervals. If they were sampled regularly, then the standard mixing trick would work:<br />\[<br />m_{n+1} = \mu m_n + (1-\mu) x_{n+1}<br />\]<br />where $m$ is our current estimate of the mean, $x_n$ is the $n$-th sample and $\mu$ determines how much history to use.<br /><br />With unequal sample times, things become a bit more complicated. If we get lots of measurements all at once, we want to give them nearly equal weight but if we have a long gap, we want to weight the very old samples much less.<br /><br />In fact, we want to weight old samples according to how old they are with exponentially decreasing weight. If we sample values $\left \lbrace x_1 \ldots x_n \right \rbrace$ at times $t_1 \ldots t_n$ then we want the weighted mean defined by<br />\[<br />m_n = {\sum_{i=1}^n x_i e^{-(t_n - t_i)/\alpha} \over \sum_{i=1}^n e^{-(t_n - t_i)/\alpha} }<br />\]<br />Here $\alpha$ plays the same role as $\mu$ did before, but on a different scale. If the evenly sampled data comes at time intervals $\Delta t$ then $\mu = e^{\Delta t / \alpha}$.<br /><br />Happily, there is a very simple recurrence relationship that allows us to keep only two intermediate values while computing the value of $m_1 \ldots m_n$ in an entirely on-line fashion as the $x_i$ values arrive.<br /><br />To see this, define<br /><br />\begin{eqnarray*}<br />\pi_n &=& e^{-(t_{n+1}-t_n)/\alpha} \\<br />w_{n+1} &=&<br />\sum_{i=1}^{n+1} e^{-(t_{n+1} - t_i)/\alpha} =<br />1+e^{-(t_{n+1}-t_n)/\alpha} \sum_{i=1}^{n} e^{-(t_{n} - t_i)/\alpha} \\<br />& =& 1 + \pi w_n\\<br />s_{n+1} &=&<br />\sum_{i=1}^{n+1} x_i e^{-(t_{n+1} - t_i)/\alpha} =<br />x_{n+1}+e^{-(t_{n+1}-t_n)/\alpha} \sum_{i=1}^{n} x_i e^{-(t_{n} - t_i)/\alpha} \\<br />&=& x_{n+1} + \pi_n s_n<br />\end{eqnarray*}<br /><br /><br />Then note that<br />\[<br />m_{n+1} = {s_{n+1} \over w_{n+1}}<br />\]<br /><br />This leads naturally to a procedure that has state consisting of $t, w, m$ which are updated with using new values of $t_n, x_n$ according to<br />\begin{eqnarray*}<br />\pi &=& e^{-(t_{n}-t)/\alpha} \\<br />w &=& 1 + \pi w \\<br />s &=& x_n + \pi s \\<br />m &=& {s \over w} \\<br />t &=& t_{n}<br />\end{eqnarray*}<br />Isn't that a kick!<br /><br />To do this right, however, we need a test. Here are some data vectors computed for $\alpha=5$:<br /><br /><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"></span><br /><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"> t x pi w s m</span><br /><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">1 11.35718 1.5992071 1.0000000 1.000000 1.5992071 1.5992071</span><br /><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">2 21.54637 -1.3577032 0.1303100 1.130310 -1.1493105 -1.0168100</span><br /><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">3 28.91061 -0.3405638 0.2292718 1.259148 -0.6040683 -0.4797436</span><br /><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">4 33.03586 0.7048632 0.4382129 1.551775 0.4401527 0.2836447</span><br /><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;">5 39.57767 0.3020558 0.2702621 1.419386 0.4210124 0.2966159</span><br /><div><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"><br /></span></div>Writing a proper unit test is an exercise left to the reader. (but I will be happy to help)Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com7tag:blogger.com,1999:blog-4453805095942812863.post-32446924023639601022010-11-04T19:59:00.000-07:002010-11-04T19:59:55.430-07:00Recent lecture on Mahout for SDForumI gave a lecture last night on recent additions to Mahout.<br /><br />The slides are here: <a href="http://www.slideshare.net/tdunning/sdforum-11042010">http://www.slideshare.net/tdunning/sdforum-11042010</a><br /><br />I can amplify this post with answers to any questions that anybody puts in the comments.Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com0tag:blogger.com,1999:blog-4453805095942812863.post-49182566844644383242010-10-31T12:23:00.000-07:002010-10-31T12:23:29.297-07:00Mahout 0.4 released!Go to the <a href="http://mahout.apache.org/">Apache Mahout</a> site for more info. Here is the official announcement:<br /><br /><blockquote>We are pleased to announce release 0.4 of Mahout. Virtually every corner of the project has changed, and significantly, since 0.3. Developers are invited to use and depend on version 0.4 even as yet more change is to be expected before the next release. Highlights include:</blockquote><blockquote>- Model refactoring and CLI changes to improve integration and consistency</blockquote><blockquote>- New ClusterEvaluator and CDbwClusterEvaluator offer new ways to evaluate clustering effectiveness</blockquote><blockquote>- New Spectral Clustering and MinHash Clustering (still experimental)</blockquote><blockquote>- New VectorModelClassifier allows any set of clusters to be used for classification</blockquote><blockquote>- Map/Reduce job to compute the pairwise similarities of the rows of a matrix using a customizable similarity measure</blockquote><blockquote>- Map/Reduce job to compute the item-item-similarities for item-based collaborative filtering</blockquote><blockquote>- RecommenderJob has been evolved to a fully distributed item-based recommender</blockquote><blockquote>- Distributed Lanczos SVD implementation</blockquote><blockquote>- More support for distributed operations on very large matrices</blockquote><blockquote>- Easier access to Mahout operations via the command line</blockquote><blockquote>- New HMM based sequence classification from GSoC (currently as sequential version only and still experimental)</blockquote><blockquote>- Sequential logistic regression training framework</blockquote><blockquote>- New SGD classifier</blockquote><blockquote>- Experimental new type of NB classifier, and feature reduction options for existing one</blockquote><blockquote>- New vector encoding framework for high speed vectorization without a pre-built dictionary</blockquote><blockquote>- Additional elements of supervised model evaluation framework</blockquote><blockquote>- Promoted several pieces of old Colt framework to tested status (QR decomposition, in particular)</blockquote><blockquote>- Can now save random forests and use it to classify new data</blockquote><blockquote>- Many, many small fixes, improvements, refactorings and cleanup</blockquote><blockquote>Details on what's included can be found in the release <a href="https://issues.apache.org/jira/secure/ReleaseNote.jspa?projectId=12310751&styleName=Html&version=12314281">notes</a>.</blockquote><blockquote>Downloads are available from the <a href="http://www.apache.org/dyn/closer.cgi/lucene/mahout/">Apache Mirrors</a></blockquote>Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com0tag:blogger.com,1999:blog-4453805095942812863.post-14268934858636385012010-10-25T11:53:00.000-07:002010-10-25T11:53:31.477-07:00New Mahout release comingThe vote has started for the 0.4 Mahout release. Lots of new stuff, but the part that I am excited about is a fairly comprehensive implementation of logistic regression suitable for large scale training and high speed classification, but there is a whole lot more.<br /><br />With the 0.4 release, Mahout is moving along strongly towards the fabled 1.0 release. At that point, we will start paying lots of attention to backwards compatibility. That will be good, but the current wild and wooly policy is pretty handy if you have something in mind that Mahout really, really needs because we can get new things in pretty readily right now.<br /><br />See http://mahout.apache.org/ for the release when it arrives and watch my twitter feed for an announcement.Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com0tag:blogger.com,1999:blog-4453805095942812863.post-38786080037447966312010-10-13T07:09:00.000-07:002010-10-13T10:16:05.553-07:00Why is the sum of two uniform randoms not uniform?Lance Norkog asks on the Mahout mailing list why adding two uniformly distributed random variables gives a pyramidal distributed value. I would normally answer on the mailing list, but here I can use lovely math notation. As I mentioned on-list, this is a very basic result that is related to the law of large numbers.<br /><br />If we were to draw a picture of the joint distribution of these variables \(x\) and \(y\), we would get something that is 1 inside the \([0,1] \times [0,1]\) square and 0 outside that region.<br /><br /><a href="http://1.bp.blogspot.com/_c4sz5uEKsbI/TLXaSGfWt6I/AAAAAAAAAFk/96oUNz0i1OI/s1600/sum-uniforms.png" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" src="http://1.bp.blogspot.com/_c4sz5uEKsbI/TLXaSGfWt6I/AAAAAAAAAFk/96oUNz0i1OI/s1600/sum-uniforms.png" /></a>For a given value \(\alpha\) of the sum \(x + y\), there is a diagonal line segment where \(x+y=\alpha\) and \(x\) and \(y\) are in the square. Where \(z \le 0\) or \(z \ge 2\) that intersection vanishes and for \(0 \lt z \lt 2\), that intersection varies in length. The probability of the sum having some particular value z is proportional to the length of that intersection. As you can imagine, the intersection varies in size linearly and it reaches a maximum where z = 1.<br /><br />For the sum of three random variables, the situation is more complex to reason about geometrically because we need to worry about the intersection of a plane and a cube. For more variables, the geometry is not worth the trouble.<br /><br />If we tackle the problem a bit more rigorously, then the easiest way to approach the problem is to compute the cumulative distribution of various values of sums. That leads to a convolution integral over the density functions involved. Since the densities are all 1, the integration limits are the key to the value and those limits have to broken down into cases. Actually doing these integrals is a pretty rare activity since the limit is approximated so well after just a few iterations.<br /><br />Just how quickly that convergence happens is can be seen by looking at the empirical distribution of the sum of three uniform deviates. I used something very like this R code to produce the graph below:<br /><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"><br /></span><br /><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"><span class="Apple-style-span" style="font-size: small;"> <b> hist(runif(100000)+runif(100000)+runif(100000), </b></span></span><br /><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"><span class="Apple-style-span" style="font-size: small;"><b> breaks=50, prob=T)</b></span></span><br /><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"><span class="Apple-style-span" style="font-size: small;"><b> lines(seq(-1,4,by=0.01), dnorm(seq(-1,4,by=0.01), </b></span></span><br /><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"><span class="Apple-style-span" style="font-size: small;"><b> 1.5, sqrt(1/4)))</b></span></span><br /><span class="Apple-style-span" style="font-family: 'Courier New', Courier, monospace;"><br /></span><br />In this graph, the red curve is the normal distribution with the same mean and standard deviation. As you can see, the peak is a tiny bit too high and the tails are just a <i>skoshi</i> too long for the normal to be a good description of the samples of the sum. <div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/_c4sz5uEKsbI/TLXbGbqVHiI/AAAAAAAAAFo/h3Ono2hqh5Q/s1600/sum_3.png" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" height="284" src="http://3.bp.blogspot.com/_c4sz5uEKsbI/TLXbGbqVHiI/AAAAAAAAAFo/h3Ono2hqh5Q/s320/sum_3.png" width="320" /></a></div><div>This is, however, just the sum of three random values. If you sum more values, the convergence to the normal distribution is very strong and by the time you are adding six uniform random values together, the difference between the distributions is no longer visible in a graph like this and can only be detected numerically using lots of data and clever things like a Kolmogorov-Smirnov test.<br /><br /><div>The moral here is that there isn't much way to avoid this regression to the normal distribution and distorting the data to avoid it is probably pointless.</div><div><br /></div><div>But if you are like me, being just a little more normal always made it easier to get along.</div></div>Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com3tag:blogger.com,1999:blog-4453805095942812863.post-69531836282503350732010-08-11T11:48:00.000-07:002010-08-11T11:48:03.848-07:00Sadder thingsI spent this last weekend holding a hand that grew cold in the early hours of Sunday morning. <br /><br />That hand helped me through much of my life. No <a href="http://www.legacy.com/obituaries/lvrj/obituary.aspx?n=hal-dunning&pid=144607613">longer</a>. At least not in the flesh.<br /><br />Nobody who reads this blog is likely to have known my father and given how little he talked about things he had done, few who knew him would know much of the many things he did. He lived a long life and a full one. Along the way he saw things few will ever see. <br /><br />In his prime, he was simply extraordinary. He could see and he could hear better than anyone I have ever known. That could be torture, as it was the time when a cat walking in the next room woke him from a deep sleep but it was what let him fly the way he did. And fly he did in planes large and small. He checked out Gary Powers in the U-2, flew P-47's and P-38 in combat and flew with me in small aircraft. We fished and walked and camped across the western US and we lived in many places. <br /><br />He didn't show off his mental abilities, but there too he could do things few others could match. He passed a graduate reading exam in French without ever studying any romance language. I saw him on several occasions understand spoken German also without having studied the language. He spoke of the shape and rate of physical processes in ways that only somebody with innate ability in math could possibly do. <br /><br />These faculties declined in age as they must with all of us, but even thus dimmed his candle burned bright.<br /><br />But it did not last. I saw it sputter and fade. Then, between one instant and the next, it was gone.Ted Dunning ... apparently Bayesianhttp://www.blogger.com/profile/02498665124454933570noreply@blogger.com3