On Bicycles, and.... what else is there?

The Golden State Warriors lost to the Cleveland Cavaliers in basketball last night. They're now down 2-1 in the series. The first team to win 4 is the champion. I don't care about basketball but there's one thing I like about the game and that's that scores generally increase relatively at random (I hope -- I hope the near miraculous comeback yesterday from a 17-point deficit after three quarters wasn't programmed), and that games are won to some degree seemingly at random. I like random. As an aside, I do find it remarkable in basketball how often a team with a big deficit claws its way back only to lose in the end by a small margin. I'd like to see a statistical analysis of this. A huge amount of money is at stake for games not being a total blow-out. I do wonder at this. Basketball has long seemed to me to be more about the show and less about a fair contest. And that makes it very difficult for me to care about the result. But this is an aside. As...

I had to do this problem for some C++ code I'm writing for work and so decided to put the result here. I have a scalar function defined over an irregular 2-dimensional mesh of triangular elements and I need to integrate the function over the surface, for example to find the average value. All I have are the coordinates of the elements. triangular mesh example, from University of Vermont Consider an element with points p 0 , p 1 , and p 2 . I can define side vectors r 0 ≡ p 1 - p 0 , and r 1 = p 2 - p 1 . Then I can determine the area of the triangle a = | r 0 × r 1 |. The integral of the function over the area is the multiplication of the average value of the function times the area. I have the area so I need the average value. This is of course ill-defined because I don't know how the function varies between the three points. But I can guess that it varies linearly. A linear function in two dimensions is defined by three degrees of freedom, for example an i...

The baccarat problem got me thinking about the random walk problem, because each of the three times in my game I reached the $200 betting limit my revenue became a random walk. At this point the problem was: which was going to come first, was I going to get into the black on $200 bets or was I going to burn through all of my cash? This is a random walk: my revenue bounces back and forth, an approximate 50% chance of each, and the game ends when I reach either of two targets: $200 above where I started (between $5000 and $6000) or $0. Many years ago, I can't remember when, I encountered a problem in random walk probabilities which was: suppose a robot, starting at x = 0, steps either in the +1 direction or the -1 direction, at random, for infinite time. What is the probability he never returns to 0 after his first step? Wow -- this is a heady problem. It always seemed to me the probability was zero: surely in infinite time he has to eventually reach all points of finite x. ...

I was watching a tutorial for some modeling software when I saw a reference to Lagrange multipliers. Lagrange multipliers... it range a bell from my distant past but I didn't recall what it was. Fortunately Google is my friend and I pulled up a video from MIT on a "recitation" -- what MIT calls sessions with teaching assistants which occur between lectures taught by professors designed to provide practice, review, and supplemental material. The problem posted on the video was the following: suppose I want to optimize (maximize or minimize) a function f(x, y, z) = x 2 + x + 2 y 2 + 3 z 2 , where the solution is constrained on the unit sphere: g(x, y, z) = x 2 + y 2 + z 2 - 1 = 0 Lagrange multipliers are based on the assumption that for points g(x, y, z) = 0 if the function f(x, y, z) is optimal along contours in g(x, y, z), then the gradient of f(x, y, z) and the gradient of g(x, y, z) must differ only by a scale factor λ, assuming the derivatives of both functi...

Last time I showed a comparison of Powertap to Vector power with a least-square line through the points. I noted this was "hard to interpret". Here's why. A traditional least-square fit begins by assuming that the x-values are perfectly known, the y-values are uncertain. It additionally assumes that the y-values are normally distributed, with a Gaussian probability distribution, with the error distribution the same for all points, or with variances inversely proportional to the weighting factors for a weighted least-square fit. This seems like a lot of technicality, but it introduces subtle biases in the result when comparing two values whcih contribute relatively equally to the error. For example, consider my power comparison for the Garmin Vector to the Powertap: The slope was 0.97. Therefore, if I flip the axis the slope should be 1.03%. If the Vector is 3% lower than Powertap than the Powertap is 3% more than Vector, on average. But that's not what I ...

Two posts back I proposed an algorithm for calculating maximal power data with a variable expiration date. For each duration, instead of retaining a single power point, I retain a list of power points where each retained point is the highest power for points of that age or newer. So going back in time for a given duration, the points are increasing power order as well as increasing age order. If my most recent effort was my highest power effort, then I have only one point in the list. In the opposite extreme, if my power decreases every activity going forward in time, I would have every activity in the list. But obviously that would be exceptional. Far more commonly powers at a given duration are sometimes higher, sometimes lower, so the number of retained points will be between the two extremes. Uploading activities to Strava on a near daily basis, it is not unusual for a user to have 1000 activities or more. So ideally this list needs to retain only a small fraction of the to...

I wanted to test the robustness of the fitting algorithm a bit more. To review, my model is the following, based on a Veloclinic model: P = P 1 ( τ 1 / t ) ( 1 - exp[-t / τ 1 ] ) + P 2 / ( 1 + t / α 2 τ 2 ) α 2 The two parameter options are: a 5-parameter model, where I fit P 1 , τ 1 , P 2 , τ 2 , and α 2 independently. a 4-parameter model, where I fit P 1 , τ 1 , P 2 , and τ 2 independently, but assume a value for α 2 . For this post, the assumed value of α 2 in the 4-parameter model was 0.5. I needed test data against to fit the model. In the first test, I used an older form of my model, where I set the coefficients the same for every test set. To generate the data I assumed a series of rides, each of constant power, each rider of a randomized duration with the durations chosen to span the data range. Each ride was of equal "quality": the power for the full duration of the "ride" was the power from the model. Fo...

introduction I spent 2 weeks analyzing fitting power-duration curves, inspired by a Veloclinic post on an analytic form for an anergy source with both rate and capacity limits (I don't know why I found this enlightening, as you get the same equation from simple capacitor-resistor networks in electrical engineering, but I never made that bridge). Anyway, Veloclinic has continued to move ahead with feedback models, but he remains committed to a model where the aerobic energy system has fixed capacity. I have stuck with more heuristic approaches where the aerobic system has time-dependent capacity. I still have a limited anaerobic capacity, but that's because attempts to replenish the anaerobic system over time would become hopelessly entangled with the aerobic power component. So simplify. I used two components with some success: an aerobic component and an anaerobic component. This is all that the CP model uses, and it oversimplifies each. Yet despite this it is the mos...

I decided to give one more try at modeling the maximal power curve in the neuromuscular power region, using the low-ceiling approach I described last time, sort of a post-hoc adjustment to calculated power which assures that what was originally a non-increasing function would remain non-increasing. First I calculate an initial power as before: P init = P 1 (τ 1 / t) (1 - exp[-t / τ 1 ]) + P 2 / (1 + t / α 2 τ 2 ) α 2 I like this model because each parameters has an easy-to-interpret meaning. It also seems to fit fairly nicely using my iteratively weighted nonlinear least square fitting scheme, where I use weights of approximately 1 for points below the envelope curve and 10 thousand for points above the envelope curve. On weighting I say "approximately 1" because I use 1 if the jump to surrounding points is 2.5%, but if the jump is more than 2.5% I scale the weight by ln( t n+1 / t n-1 ) / 2 ln 1.025. Note I don't apply this to the points falling abov...

I've noted that the model persistently predicts, for me, higher sprint powers than I have recorded. Is this difference significant? Perhaps I simply tend to underperform relative to capacity for very short sprints. Or perhaps it is an anomaly associated with the powertap power measurement (these data were mostly taken on the yellow Cervo computer, by the way). But assuming it is significant, it would be preferable to try and model it. The key assumption is that each component is able to deliver its full capacity all the way to zero time. This is an obvious assumption for any given component, but since the model simply adds components, the implicit assumption is that I am able to deliver my maximum contribution from all components simultaneously. So if I am able to deliver my maximum component-1 power, which I call anaerobic, I am simultaneously able to deliver my maximum component 2 power, which I'll call aerobic. Since aerobic power can be sustained for a long time (n...

Last time I proposed the following functional form for the aerobic component: f = 1 / [ 1 + (t / τ)] α For times much less than the threshold τ this saturates to 1. For times much longer, than τ it has a power-law dependence with exponent α. Differentiating this with respect to time yields the following: ∂ f / ∂ t = -(α / τ) / [1 + (t / τ)] α + 1 In the limit of small t this has the value -α / τ, so no singularity. This would add an additional fitting parameter to the model. A reasonable initial guess for α might be 0.5. One issue here is that the time derivative, to first order, is is proportional to α / τ, while the value for small times depends on neither to zero order. So if the fitting algorithm is working separately on α and τ, it will have a difficult time separating the contribution of α from the contribution of τ. It would be tempted to increase or decrease both tog...

Previously I fit a 2-component "modified Veloclinic" model to some of my old power-duration data. This worked fairly well, but sprinting power was almost universally overestimated. Experimentally there was a sharper transition from time-dependent power to peak sprint power. The thing is I don't see this in Veloclinic's fits. Nor do I see it in the Damien's Gauser's enhanced model in Golden Cheetah 3.1 development version, which is curious. So what's the problem? Well, one difference is I embraced the strictly heuristic modification of the aerobic term with an α term: P aerobic = P 2 ( τ 2 / t) α 2 ( 1 - exp[ -( t / τ 2 ) α 2 ] ) Where I observed that for this term α = 1/2 worked nicely. So I decided to go with that as a constant, and not fit it to the models. I got better fits to the long-time data than I see with Veloclinic's model, which has α = 1 for all components. The rationale is that aerobic power isn't fixed capacity, but h...

Previously I used a damping factor δ for the nonlinear fitting of the modified Veloclinic power-duration model. This resulted in me fitting 28 out of 29 test cases when I set δ = 0.25. But δ = 0.25 is too small: it means at best 25% of the progress toward the solution is covered with each iteration. That's inefficient when homing in on the result. The real issue is when the solver isn't homing in on the result: when the solution is still substantially off and the solver it attempting to make a big leap toward the solution. The big leap could easily overshoot the desired solution, or be in a slightly wrong direction in the hyperdimensional parameter space. In such instances, it's better to take smaller steps toward the solution, to home in on the the "zone of quadratic convergence" where the targeting becomes easier. For that I introduced a dynamic weighting scheme multiplying the primary damping factor by a factor dependent on the step size of the parameters. ...

Wanting to do some matrix math, I installed the Math::Matrix module from CPAN . There's some nice instructions on installing CPAN modules on about.com . A common tool people use for matrix manipulation is MATLAB . I'm deeply embarrassed I have virtually no skills in MATLAB. However, this module is pretty cool. What makes it so nice it is uses overloading of the basic mathematical operators. So that makes it easy to add, subtract, and multiply matrices, as well as a string conversion when it's accessed in a scalar context. Here's an example: use Math::Matrix; use strict; my $a = new Math::Matrix ([1, 0, -1], [-1, 1, 0], [0, -1, 1]); warn("a =\n$a\n"); my $at = $a->transpose; warn("at =\n$at\n"); warn("at * a =\n", $at * $a, "\n"); The result: a = 1.00000 0.00000 -1.00000 -1.00000 1.00000 0.00000 0.00000 -1.00000 1.00000 at = 1.00000 -1.00000 0.00000 0.00000 1.00000 -1.0000...

Two posts ago I described fitting the linear CP model with an iteratively weighted least-square fit to approximate an envelope fit. The weights were either 1/t 2 for points falling under the CP curve or 10k/t 2 for points falling above the CP curve. This did a decent job in the example I showed of matching the fit using a 2-point method used in Golden Cheetah. But of greater interest is using this method for my modified Veloclinic model. That model should do much better at fitting the full time spectrum of the power-duration data. The challenge there, however, is the model is nonlinear. Nonlinear least squares fitting is described by Wolfram . In my last posts, I added weighting to the method described on the Wolfram page. Nonlinear least square fitting is a matter of navigating a hyperdimensional space, looking for the point where all error terms go to zero (where the model perfectly fits the data). There is no such point, so you search for it until improvement slows below s...

Note: this uses MathML code produced by OpenOffice. It doesn't seem to work on Internet Explorer. It's also been rejected by Chrome . But it works on Safari, Chromium, and Firefox. Nonlinear least-square fits are done with a Jacobian Matrix, which describes the linearized dependencies of the model, evaluated at each of the points, on each of the model parameters.  In this case I show two parameters, τ 1 and τ 2 .  I like writing these things out, rather than using index notation, because index notation is a bit abstract.  In the following, f is a function of time t representing the model to be fit. A = ∂ f ∂ τ 1 ∣ t 1 ∂ f ∂ τ 2 ∣ t 1 ∂ f ∂ τ 1 ∣ t 2 ∂ f ∂ τ 2 ∣ t 2 ⋮ ⋮ ∂ f ∂ τ 1 ∣ t N − 1 ∂ f ∂ τ 2 ∣ t N − 1 ∂ f ∂ τ 1 ∣...