//==================================================================
// Carl Schmertmann
// created 31 Jul 2004
// altered 02 Mar 2005
//
// JavaScipt source code for NLS estimation of Quadratic Spline
// fertility schedules
//==================================================================
//-- GLOBAL VARIABLES
var CR = unescape("%0D");
var LF = unescape("%0A");
var TB = unescape("%09");
var NL = CR + LF;
var xval = new Array()
var fval = new Array()
var nval = new Array()
// xx and yy are for plotting
var xx = new Array()
var yy = new Array()
var xval_text = new Array();
var fval_text = new Array();
var DataReport = ""
var IterReport = ""
var FitNReport = ""
var Fit1Report = ""
var nint = 999 // # of age intervals in input data supplied by user
var DataGrabbed = false
var MaxIter = 50
var R = 0; var alpha = 0; var P = 0 ; var H = 0 ; var beta = 0;
var RDefault = ".200"
var alphaDefault = "14"
var PDefault = "23"
var HDefault = "36"
//-- example data: period 5fx values from Iran 2002
var ExampleData = "15" + TB + " 19.3" + NL +
"20" + TB + " 86.8" + NL +
"25" + TB + "136.4" + NL +
"30" + TB + "100.5" + NL +
"35" + TB + " 42.5" + NL +
"40" + TB + " 15.1" + NL +
"45" + TB + " 2.1" ;
//-- example data: period 1fx values from Sweden 1996
var ExampleDataSingleYear = "15" + TB + "0.000547" + NL + "16" + TB + "0.002131" + NL + "17" + TB + "0.005386" + NL +
"18" + TB + "0.009868" + NL + "19" + TB + "0.020513" + NL + "20" + TB + "0.032461" + NL + "21" + TB + "0.044743" + NL +
"22" + TB + "0.056310" + NL + "23" + TB + "0.071642" + NL + "24" + TB + "0.083872" + NL + "25" + TB + "0.103400" + NL +
"26" + TB + "0.112203" + NL + "27" + TB + "0.119024" + NL + "28" + TB + "0.120221" + NL + "29" + TB + "0.121122" + NL +
"30" + TB + "0.113856" + NL + "31" + TB + "0.105922" + NL + "32" + TB + "0.092814" + NL + "33" + TB + "0.081299" + NL +
"34" + TB + "0.069526" + NL + "35" + TB + "0.058160" + NL + "36" + TB + "0.049225" + NL + "37" + TB + "0.037610" + NL +
"38" + TB + "0.028761" + NL + "39" + TB + "0.021161" + NL + "40" + TB + "0.014638" + NL + "41" + TB + "0.009330" + NL +
"42" + TB + "0.006259" + NL + "43" + TB + "0.003281" + NL + "44" + TB + "0.001776" + NL + "45" + TB + "0.000881" + NL +
"46" + TB + "0.000325" + NL + "47" + TB + "0.000156" + NL + "48" + TB + "0.000031" + NL + "49" + TB + "0.000000" ;
//-- 2 to 1 dim index translation: j and k are row and column in an NRxNC matrix, result is matching index in an (NRxNC)x1 vector
function ix(j,k,NR,NC) { return NC*j +k }
//-- pretty output: fixed the number of decimal places D in a string of width W<12
function Pretty(x,W,D) {
st = " " + x.toFixed(D)
st = st.substr(st.length-W , W)
return st
}
//-- pretty string output: right-justified string of width W<12
function PrettyStr(s,W) {
st = " " + s
st = st.substr(st.length-W , W)
return st
}
//-- Reset the form to default values
function RestoreDefault(form) {
form.R_text.value = RDefault
form.alpha_text.value = alphaDefault
form.P_text.value = PDefault
form.H_text.value = HDefault
form.data.value = ""
form.ReportType.value = "Data"
form.report.value = ""
DataReport = ""
IterReport = ""
FitNReport = ""
Fit1Report = ""
PlotReset()
document.plotapp.SetHeader("QS Fit")
DataGrabbed = false
} ;
//-- Load the example data (Iran 2002 period 5fx values)
function LoadExampleData(form,h) {
if ( h==5 ) { form.data.value = ExampleData } else { form.data.value = ExampleDataSingleYear }
PlotReset()
DataGrabbed = false
} ;
//-- Get the (R,alpha,P,H) values from the screen text and put into numeric variables
function GetParms(form) {
R = eval( form.R_text.value)
alpha = eval( form.alpha_text.value)
P = eval( form.P_text.value)
H = eval( form.H_text.value)
}
//-- Get the (x,nfx) values from the left-hand text area of the form and put into numeric variables xval and fval
//-- this function is adapted from TW Shattuck, http://www.colby.edu/chemistry/PChem/scripts/lsfitpl.html
function GetData(form) {
var line = new Array() ;
var instr = new Array(4) ;
var rews = /[\s\t,]+/g ; // regular expression matching white space
var resws = /^\s+/ ; // matches only white space at the beginning of the line
var lstart = 1 ; var lend = 4 ;
var maxf = -999.0;
var maxix = 0;
var ix = 0;
line = form.data.value.split("\n") ; // splits text area into line[0], line[1], ...
if ( line.length <= 1 ) line = form.data.value.split("\r") ;
nint = parseInt(line.length) ;
for (i=0 ; i maxf) {
maxf = fval[i]
maxix = i
}
if ( isNaN(xval[i]) ) { alert("Input error x="+instr[0]) ; return false }
if ( isNaN(fval[i]) ) { alert("Input error y="+instr[1]) ; return false }
} ;
nint = i ; // important variable: this is the number of intervals for which we're fitting rates
if (!DataGrabbed) {
R = maxf;
P = (xval[maxix] + xval[maxix+1])/2
form.R_text.value = R.toFixed(3) // scale initial R value according to input data
form.P_text.value = P.toFixed(2) // set initial P value according to input data
ix = maxix + 1
while ((fval[ix] > maxf / 2) && (ix < nint-1)) { ix += 1 }
if (ix < nint-2) { H = (xval[ix] + xval[ix+1])/2 } else { H = xval[ix] }
form.H_text.value = H.toFixed(2) // set initial H value according to input data
}
document.plotapp.SetXYRange(10,50,0, 1.33*maxf) // set the scale for f(x)
// cycle through one more time and calculate interval widths nval[0]...nval[nint-1]
for (i=0 ; i 50) { beta = Z1} ; if (Z2 < 50) { beta = Z2}
knot[0] = alpha
knot[1] = (1-W)*alpha + W*P
knot[2] = P
knot[3] = (P+H)/2
knot[4] = (H+beta)/2
theta[0] = 1 / (W * (P-alpha) * (P-alpha)) ; theta[1] = -theta[0] / (1-W)
ZA = (H-knot[2]) * (H-knot[2]) ; ZB = (H-knot[3]) * (H-knot[3])
ZC = (beta-knot[2]) * (beta-knot[2]) ; ZD = (beta-knot[3]) * (beta-knot[3]) ; ZE = (beta-knot[4]) * (beta-knot[4])
ZF = 2 * (beta-knot[2]) ; ZG = 2 *(beta-knot[3]) ; ZH = 2 *(beta-knot[4])
Z1 = 0.5 - (theta[0] * (H-alpha)*(H-alpha) + theta[1] * (H-knot[1])*(H-knot[1]))
Z2 = 0 - (theta[0] * (beta-alpha)*(beta-alpha) + theta[1] * (beta-knot[1])*(beta-knot[1]))
Z3 = 0 - (2*theta[0] * (beta-alpha) + 2*theta[1] * (beta-knot[1]))
DENOM = ZA*ZD*ZH - ZA*ZE*ZG - ZC*ZB*ZH + ZF*ZB*ZE
theta[2] = (Z1*(ZD*ZH-ZE*ZG) - Z2*(ZB*ZH) + Z3*(ZB*ZE) ) / DENOM
theta[3] = (Z1*(ZE*ZF-ZC*ZH) + Z2*(ZA*ZH) - Z3*(ZA*ZE) ) / DENOM
theta[4] = (Z1*(ZC*ZG-ZD*ZF) + Z2*(ZB*ZF-ZA*ZG) + Z3*(ZA*ZD-ZB*ZC) ) / DENOM
for (i=0; i < 5; i++) { ResultVec[i] = theta[i] ; ResultVec[5+i] = knot[i] }
ResultVec[10] = beta
return ResultVec
}
//-- calculate predicted-observed fertility rates for the age intervals in the input data, as a function of (R,alpha,P,H)
//-- this uses the exact analytical expressions for CUMULATIVE fertility, CumulF(x) = R/3 * sum { theta[k] * [max(0,x-t[k])]^3 }
function Resid(R,alpha,P,H) {
var epsilon = Array();
var theta = Array();
var knot = Array();
var beta = 50.0
//alert("in Resid, R = " + R + " alpha= " + alpha + " P= " + P + " H= " + H)
sp = SplineParams(R,alpha,P,H)
for (i=0; i < 5; i++) { theta[i] = sp[i] ; knot[i] = sp[5+i] }
beta = sp[10]
for (i=0; i < nint ; i++) {
Fa = 0
Fb = 0
a = Math.min(Math.max(0, xval[i]), beta)
b = Math.min(Math.max(0, xval[i]+nval[i]), beta)
for (k=0 ; k < 5 ; k++) {
za = Math.max( 0 , a - knot[k])
zb = Math.max( 0 , b - knot[k])
Fa += R/3 * theta[k] * za * za * za
Fb += R/3 * theta[k] * zb * zb * zb
}
epsilon[i] = (Fb-Fa)/(b-a) - fval[i]
}
return epsilon
}
//-- calculate predicted fertility at exact ages in the xx array -- for plotting
function Predict(xx,R,alpha,P,H) {
var epsilon = Array();
var theta = Array();
var knot = Array();
var beta = 50.0
sp = SplineParams(R,alpha,P,H)
for (i=0; i < 5; i++) { theta[i] = sp[i] ; knot[i] = sp[5+i] }
beta = sp[10]
for (i=0; i < xx.length ; i++) {
fx = 0
if (xx[i] < beta) {
for (k=0 ; k < 5 ; k++) {
z = Math.max( 0 , xx[i] - knot[k])
fx += R * theta[k] * z * z
}
}
yy[i] = fx
}
return yy
}
//-- calculate TFR
function TFR(R,alpha,P,H) {
var theta = Array();
var knot = Array();
var beta = 50.0
sp = SplineParams(R,alpha,P,H)
for (i=0; i < 5; i++) { theta[i] = sp[i] ; knot[i] = sp[5+i] }
beta = sp[10]
Fb = 0
b = beta
for (k=0 ; k < 5 ; k++) {
zb = Math.max( 0 , b - knot[k])
Fb += R/3 * theta[k] * zb * zb * zb
}
return Fb
}
//-- calculate predicted 1fx values at ages 10...49
//-- the first 40 elements of the output vector contain rates for one-year intervals [1f10...1f49]
//-- the next 40 elements contain rates at exact ages [f(10)...f(49)]
function Fit1fx(R,alpha,P,H) {
var fsingle = Array();
var theta = Array();
var knot = Array();
var beta = 50.0
sp = SplineParams(R,alpha,P,H)
for (i=0; i < 5; i++) { theta[i] = sp[i] ; knot[i] = sp[5+i] }
beta = sp[10]
for (x=10; x < 50 ; x++) {
Fa = 0
Fb = 0
exactF = 0
a = Math.min(Math.max(0, x), beta)
b = Math.min(Math.max(0, x+1), beta)
for (k=0 ; k < 5 ; k++) {
za = Math.max( 0 , a - knot[k])
zb = Math.max( 0 , b - knot[k])
Fa += R/3 * theta[k] * za * za * za
Fb += R/3 * theta[k] * zb * zb * zb
exactF += R * theta[k] * za * za
}
if (x > beta) { exactF = 0 }
fsingle[x-10] = (Fb-Fa)
fsingle[x-10+40] = exactF
}
return fsingle
}
//-- Solve a matrix equation A x = b for x. A must be 4x4.
//-- The single 4x5 argument to this function is [A | b] going in, transformed to [inv(A) | x] coming out
//-- This function is adapted from John Pezzullo, http://members.aol.com/johnp71/nonlin.html
function Msolve(Z) {
var s = 0.0
for (i=0 ; i < 4 ; i++) {
s = Z[ix(i,i,4,5)]
Z[ix(i,i,4,5)] = 1.0
for (k=0 ; k < 5; k++) { Z[ix(i,k,4,5)] = Z[ix(i,k,4,5)] / s }
for (j=0; j < 4; j++) {
if (i !=j) {
s = Z[ ix( j,i,4,5)];
Z[ ix( j,i,4,5)] = 0;
for (k=0; k < 5; k++) {
Z[ ix(j,k,4,5)] = Z[ ix(j,k,4,5)] - s * Z[ ix(i,k,4,5) ]
}
}
}
}
}
//-- PLOTTING FUNCTIONS
//-- blank out all the lines on the plot (including bar graph of data)
function PlotReset() {
document.plotapp.ResetMeter()
}
//-- add the bar graph for the current (x, f) data -- this will be line "1" on the XY plot
function PlotBars(xval,fval,nval) {
for (i=0; i < nint; i++) {
if (i < nint-1) { xnext=xval[i+1] } else { xnext = 50 }
document.plotapp.SetMeterXY(1, xval[i], 0)
document.plotapp.SetMeterXY(1, xval[i], fval[i])
document.plotapp.SetMeterXY(1, xnext, fval[i])
document.plotapp.DrawMeterXY(1, xnext, 0)
}
}
//-- add a QS line -- this will be line "2" on the XY plot
function PlotLine(xx,yy) {
document.plotapp.ResetMeter() // blank out current lines
PlotBars(xval,fval,nval) // add bars back
// redraw QS line
npts = xx.length
for (i=0; i < npts-1; i++) { document.plotapp.SetMeterXY(2, xx[i], yy[i]) }
document.plotapp.DrawMeterXY(2,xx[npts-1],yy[npts-1])
}
//-- FIT THE MODEL
function QSFit(form) {
var eps = Array()
var p = Array()
var deps = Array()
var delta = [ 0.00005 , 0.0005 , 0.0005, 0.0005 ]
var Z = Array()
var Q = Array()
var scale = [ 0, 0.01, 0.10, 0.5, 1, 2 , 5, 10 ]
var best_index = 0
var bestQ = Infinity
var ee = Array()
var SearchDirection = Array()
var stderr = Array()
var SSE = 0.0
var fsingle = Array()
function ResetParms() { p[0] = R ; p[1] = alpha ; p[2] = P ; p[3] = H }
function BadParms(R,alpha,P,H) { return (alpha > P) || (P > H) }
//-- define 'convergence' as less than 1 part per 10,000 change in objective function & all parameters
function convergence() {
Qstar = bestQ
conv = true
conv = (conv && ( Math.abs( (Qstar - SSE)/SSE) < .0001 ) )
conv = (conv && ( Math.abs( (bestR-R) / R ) < .0001 ))
conv = (conv && ( Math.abs( (bestalpha-alpha) / alpha ) < .0001 ))
conv = (conv && ( Math.abs( (bestP-P) / P ) < .0001 ))
conv = (conv && ( Math.abs( (bestH-H) / H ) < .0001 ))
return conv
}
function Deriv() {
// calculate residuals at current parameter values
eps = Resid(R,alpha,P,H)
// calculate numerical derivs of residuals[0...n] with respect to (R,alpha,P,H)
for (j=0 ; j<4 ; j++) {
ResetParms()
p[j] = p[j] + delta[j]
D = Resid(p[0], p[1], p[2] , p[3])
for (i=0 ; i= SSE) {
alert("DID NOT FIND IMPROVEMENT IN SSE -- STOPPING")
break
}
// otherwise reset parameters to best values found
else {
R = bestR
alpha = bestalpha
P = bestP
H = bestH
form.R_text.value = R.toFixed(4)
form.alpha_text.value = alpha.toFixed(2)
form.P_text.value = P.toFixed(2)
form.H_text.value = H.toFixed(2)
}
if (iter==(MaxIter-1)) { alert("EXCEEDED MAXIMUM OF " + MaxIter + " ITERATIONS -- STOPPING") }
} // for iter
} // end QSFit
//-- initialize the xx array for plots
for (age=10 ; age <= 50; age++) { xx[age-10] = age }