Log in
Register|
ElMaestro ★★★ Denmark, 2011-11-18 09:11 (4921 d 22:40 ago) (edited on 2011-11-18 09:58) Posting: # 7680 Views: 5,556 |
|
|
Hi experts, The short version: How do you derive equations such as the ones Potvin et al. gave for calculation of the residual variance from a normal linear model? The long version (you definitely want to skip this): I would like to go beyond Potvin's designs and explore a bit such as parallel, imbalance and that sorta of stuff. Set up some code in R, but it is very clear that simulating and evaluating a million trials with e.g. N1=30 takes forever on my computer and 1. I am not equipped with great patience 2. I do not have access to super-computers. A faster way of getting to the residual sigma is necessary. In my desperation I had a brilliant idea after googling a bit around. The normal linear model is y=Xb+e and the solution for the effects is b=(XtX)-1Xty and the SOS comes out simply through ete with e=Y-Xb. So I hardcoded a light-weight matrix library in C (only needed it to be able to do transpose, subtract, invert, multiply, so there is no need for the heavy stuff like GSL) to do the trick and it works perfectly. But... although it is a bit faster than R, it is just not fast at all and I can't optimise it further. So I am abandoning this idea now. Hence, I am on the prowl for simpler equations for the residual sigma. Any ideas? Anyone able to tell a person with a walnut-sized brain how the equations from Potvin can be derived and extended to more complex scenarios? Best regards, EM. PS: Any wise guy wanting to tell me I shouldn't find b through inversion but use a QR factorisation in stead? |
|
d_labes ★★★ Berlin, Germany, 2011-11-18 16:44 (4921 d 15:07 ago) @ ElMaestro Posting: # 7681 Views: 4,617 |
|
|
Dear Supreme Grand Maestro, ❝ PS: Any wise guy wanting to tell me I shouldn't find b through inversion but use a QR factorisation in stead? I must confess that I'm surely not an expert in numerical Al-Gore-Rhythms optimized for speed for the solution of a general linear model fit. And my last encounter with matrices calculus is also some long years ago. Nowadays I had only coincidences with Matratzen every night  .And I'm surely not wise even so old. But just some small observations from what R itself is doing: Type in your R console lm without brackets (!) to see in the code that the working horse of this function is the function lm.fit().?lm.fit will give you the information that it uses a QR factorisation as the solely method for solving. Type lm.fit (again without brackets) and notice that here the calculation burden is assigned to a call of compiled FORTRAN code via function .Fortran(...), namely to the LINPACK routine DQRDC as the help tells us. The rest are some native R code lines surely not important for speed.Thus I highly suspect that you don't get great speed gains if you implement the whole process in FORTRAN or C by yourself. — Regards, Detlew |
|
ElMaestro ★★★ Denmark, 2011-11-20 19:54 (4919 d 11:57 ago) @ d_labes Posting: # 7683 Views: 4,482 |
|
|
Hmmmm, I got a very cool idea this morning: X is invariant for any given N1 and N2 (let those be the number of experimental units at stage 1 and 2). So (XtX)-1Xt is invariant for any given N1 and N2. For any givev N1, I can generate all possible (XtX)-1Xt's and cache them in memory, just keeping track of the pointers. It consumes sh!tload of memory to go up to e.g. N2=1000, but hey so does Firefox or Excel or whatever. So, once we have y for N1 and an Algore Rhytm for stopping or proceeding to stage 2, we can just calculate N2 and, sim another N2 y-values and append to the first y-vector, then access (XtX)-1Xt that applies that that (N1+N2) through the pointers and multiply it onto the simulated combined y-vector. This is a quite miraculous speed-up. One million sims in three seconds and so far all solutions that I checked agree with R's lm.Now, one of the laws of nature state that when ElMaestro gets a good idea it is usually a bad idea after closer inspection. So I better get started digging into this  Wishing yawl a pleasant rest of the weekend. EM. |
Ing. Helmut Schütz