Nothing
if(!require("numDeriv"))stop("this test requires numDeriv.")
Sys.info()
#######################################################################
# Test gradient and hessian calculation in genD using data for calculating
# curvatures in Bates and Watts.
#model A p329,data set 3 (table A1.3, p269) Bates & Watts (Puromycin example)
#######################################################################
puromycin <- function(th){
x <- c(0.02,0.02,0.06,0.06,0.11,0.11,0.22,0.22,0.56,0.56,1.10,1.10)
y <- c(76,47,97,107,123,139,159,152,191,201,207,200)
( (th[1] * x)/(th[2] + x) ) - y
}
D.anal <- function(th){
# analytic derivatives. Note numerical approximation gives a very good
# estimate of these, but neither give D below exactly. The results are very
# sensitive to th, so rounding error in the reported value of th could explain
# the difference. But more likely th is correct and D has been rounded for
# publication - and the analytic D with published th seems to work best.
# th = c(212.70188549 , 0.06410027) is the nls est of th for BW published D.
x <- c(0.02,0.02,0.06,0.06,0.11,0.11,0.22,0.22,0.56,0.56,1.10,1.10)
y <- c(76,47,97,107,123,139,159,152,191,201,207,200)
cbind(x/(th[2]+x), -th[1]*x/(th[2]+x)^2, 0, -x/(th[2]+x)^2, 2*th[1]*x/(th[2]+x)^3)
}
# D matrix from p235. This may be useful for rough comparisons, but rounding
# used for publication introduces substantial errors. check D.anal1 - D.BW
D.BW <- t(matrix(c(
0.237812, -601.458, 0, -2.82773, 14303.4,
0.237812, -601.458, 0, -2.82773, 14303.4,
0.483481, -828.658, 0, -3.89590, 13354.7,
0.483481, -828.658, 0, -3.89590, 13354.7,
0.631821, -771.903, 0, -3.62907, 8867.4,
0.631821, -771.903, 0, -3.62907, 8867.4,
0.774375, -579.759, 0, -2.72571, 4081.4,
0.774375, -579.759, 0, -2.72571, 4081.4,
0.897292, -305.807, 0, -1.43774, 980.0,
0.897292, -305.807, 0, -1.43774, 980.0,
0.944936, -172.655, 0, -0.81173, 296.6,
0.944936, -172.655, 0, -0.81173, 296.6), 5,12))
cat("\nanalytic D:\n")
print( D.anal <- D.anal(c(212.7000, 0.0641)), digits=16)
cat("\n********** note the results here are better with d=0.01 ********\n")
cat("\n********** in both relative and absolute terms. ********\n")
cat("\nnumerical D:\n")
print( D.calc <- genD(puromycin,c(212.7000, 0.0641), method.args=list(d=0.01)),
digits=16)
# increasing r does not always help
#D.calc <- genD(puromycin,c(212.7000, 0.0641), r=10)#compares to 0.01 below
#D.calc <- genD(puromycin,c(212.7000, 0.0641), d=0.001)
cat("\ndiff. between analytic and numerical D:\n")
print( D.calc$D - D.anal, digits=16)
cat("\nmax. abs. diff. between analtic and numerical D:\n")
print( max(abs(D.calc$D - D.anal)), digits=16)
# These are better tests except for 0 column, so add an epsilon
cat("\nrelative diff. between numerical D and analytic D (plus epsilon):\n")
print(z <- (D.calc$D - D.anal) / (D.anal + 1e-4), digits=16)
# d=0.0001 [12,] 1.184044172787111e-04 7.451545953037876e-03
# d=0.01 [12,] 1.593395089728741e-08 2.814629092064831e-07
cat("\nmax. abs. relative diff. between analtic and numerical D:")
print( max(abs(z)), digits=16)
if(max(abs(z)) > 1e-6) stop("BW test FAILED")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.