from numpy import * #fitting world population data to modelling function f(t) = (t-1960) + x0 + x1 * pow(1.02,t) #that is, we assume growth rate of 1.02 (which was accurate for 1960s) #create matrix A with entries (i, pow(1.02,i)) D = []; for i in range(0,10): D = D + [[i, pow(1.02,i)]]; A = array(D); #population data from http://www.census.gov/ipc/www/idb/worldpop.php #fullb contains entries from 1960 through 2050 (actual up to 2009) #we'll use values from 1960 through 1969 to interpolate later values b = array([3041697768,3082701189,3138922062,3208676975,3280050222,3349287332,3419728217,3489679443,3561729216,3636566653]); fullb = array([3041697768, 3082701189, 3138922062, 3208676975, 3280050222, 3349287332, 3419728217, 3489679443, 3561729216, 3636566653, 3711961664, 3789539409, 3865804060, 3941550996, 4016055456, 4088611586, 4159762924, 4231510486, 4303133652, 4377496958, 4452547522, 4528882162, 4608681758, 4690277509, 4770468296, 4852051940, 4935873592, 5022023110, 5108859812, 5195713196, 5283687429, 5367185126, 5451671517, 5534137903, 5615310526, 5696676739, 5776856845, 5855086722, 5932090951, 6008254795, 6083550220, 6159242510, 6234361771, 6308516183, 6382433547, 6456443080, 6530984631, 6606214786, 6681112529, 6755987239, 6830586985, 6905172554, 6979312828, 7053079305, 7126724753, 7200008584, 7272849777, 7345201125, 7416852578, 7487663656, 7557514266, 7626477056, 7694592508, 7761728812, 7827803562, 7892758722, 7956667160, 8019596128, 8081508492, 8142383845, 8202205367, 8261039699, 8318955615, 8375945217, 8431992573, 8487066951, 8541209793, 8594470749, 8646831011, 8698265300, 8748743446, 8798303380, 8846986545, 8894758313, 8941579812, 8987418540, 9032303796, 9076268030, 9119291753, 9161358865, 9,202,458,484 ]); #create normal system Ap x = bp by multiplying with transpose of A Ap = dot(A.transpose(), A); bp = dot(A.transpose(), b); #estimate for parameters x0 and x1 x = linalg.solve(Ap, bp); print("(%r, %r)"%(x[0],x[1])); #compare to actual values for i in range(0,50): year = 1960 + i; prepop = int(x[0] + x[1] * pow(1.02,i)); actual = int(fullb[i]); error = 1.0 * (prepop - actual)/actual; print("Year: %r"%year+" Predict Pop: %i"%prepop+" Actual: %i"%actual+" Error: %r"%error );