Auto-constructing Fisher Information matrix

Description

Auto-constructs Fisher information matrix for nonlinear and generalized linear models as two R functions.

Usage

cfisher(ymean, yvar, ndpoints, prec = 53)

Arguments

Details

If response variables have the same constant variance, for example σ^2, then yvar must be 1.

Value

a list containing two closures:

For more details, see 'Note'.

Note

This function is applicable for models that can be written as E(Y_i)=f(x_i,β) where y_i is the ith response variable, x_i is the observation vector of the ith explanatory variables, β is the vector of parameters and f is a continuous and differentiable function with respect to β. In addition, response variables must be independent with distributions that belong to the Natural exponential family. Logistic,Poisson, Negative Binomial, Exponential, Richards, Weibull, Log-linear, Inverse Quadratic and Michaelis-Menten are examples of these models.

Consider a p-parameter model and a design ξ that contains n m-dimensional points. Then

x = (x_1, x_2, …, x_i, …, x_n),

w = (w_1, w_2,…, w_n),

β = (β_1, β_2, …, β_p),

where x_i = (x_{i1}, x_{i2}, …, x_{im}) is the ith design point.

Author(s)

Ehsan Masoudi, Majid Sarmad and Hooshang Talebi

References

Masoudi, E., Sarmad, M. and Talebi, H. 2012, An Almost General Code in R to Find Optimal Design, In Proceedings of the 1st ISM International Statistical Conference 2012, 292-297.

Examples

## Logistic dose response model
ymean <- "(1/(exp(-b2 * (x1 - b1)) + 1))"
yvar <- "(1/(exp(-b2 * (x1 - b1)) + 1)) * (1 - (1/(exp(-b2 * (x1 - b1)) + 1)))"
res <- cfisher(ymean, yvar, ndpoints = 2, prec = 54)

# res$fim is Fisher information matrix for a two-points design 
res$fim(x = c(x11 = 2, x21 = 3), w = c(w1 = .5, w2 = .5), b = c(b1 = .9, b2 = .8))

# res$fim is Fisher information matrix for a two-points design with 54 precision
res$fim.mpfr(x = c(x11 = 2, x21 = 3), w = c(w1 = .5, w2 = .5), b = c(b1 = .9, b2 = .8))

# Fisher information matrix for model:
fim<- cfisher(ymean, yvar, ndpoints = 1, prec = 54)
res$fim(x = c(x11 = 2), w = c(w1 = 1), b = c(b1 = .9, b2 = .8))

## posison with E(y) = exp(b1 + b2 * x1 + b3 * x1^2 + b4 * x2 +b5 * x2^2 + b6 * x1 * x2)
ymean <- yvar <- "exp(b1 + b2 * x1 + b3 * x1^2 + b4 * x2 +b5 * x2^2 + b6 * x1 * x2)"
fim <- cfisher(ymean, yvar, ndpoints = 6, prec = 54)

# res$fim is Fisher information matrix for a six-points design 
res$fim(x = c(1:12), w = rep(1/6, 6), b = c(1:6)) ## NAN

# res$fim.mpfr is Fisher information matrix for a six-points design with 53 precision
res$fim.mpfr(x = c(1:12), w = rep(1/6, 6), b = c(1:6))

## Linear regression with two indeoendent varibales (the design points are two-dimensional)
ymean <- "b1 + b2 * x1 + b3 * x2"
yvar = "1"
res <- cfisher(ymean, yvar, ndpoints = 3, prec = 54)
res$fim(x = c(1:6), w = c(.3, .3, .3))
res$fim.mpfr(x = c(1:6), w = c(.3, .3, .3))

## Logistic model:
ymean <- "1/(exp(-b1 - b2 * x1) + 1)"
yvar <- "(1/(exp(-b1 - b2 * x1) + 1)) * (1 - (1/(exp(-b1 - b2 * x1) + 1)))"
cfisher(ymean, yvar, ndpoints = 2, prec = 54)

## Poisson model:
ymean <- yvar <-  "exp(b1 + b2 * x1)"
cfisher(ymean, yvar, ndpoints = 2, prec = 54)

## Poisson dose response model:
ymean <- yvar <- "b1 * exp(-b2 * x1)"
cfisher(ymean, yvar, ndpoints = 2, prec = 54)

## Inverse Quadratic model:
ymean <- "(b1 * x1)/(b2 + x1 + b3 * (x1)^2)"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 3, prec = 54)
#
ymean <- "x1/(b1 + b2 * x1 + b3 * (x1)^2)"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 3, prec = 54)

## Weibull model:
ymean <- "b1 - b2 * exp(-b3 * x1^b4)"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 4, prec = 54)

## Richards model:
ymean <- "b1/(1 + b2 * exp(-b3 * x1))^b4"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 4, prec = 54)

## Michaelis-Menten model:
ymean <- "(b1 * x1)/(1 + b2 * x1)"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 2, prec = 54)
#
ymean <- "(b1 * x1)/(b2 + x1)"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 2, prec = 54)
#
ymean <- "x1/(b1 + b2 * x1)"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 2, prec = 54)

## log-linear model
ymean <- "b1 + b2 * log(x1 + b3)"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 3, prec = 54)

## Exponential model:
ymean <- "b1 + b2 * exp(x1/b3)"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 3, prec = 54)

## Emax model:
ymean  <- "b1 + (b2 * x1)/(x1 + b3)"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 3, prec = 54)

## Negative binomial model Y ~ NB(E(Y), theta) where E(Y) = b1*exp(-b2*x1):
theta = 5
ymean <- "b1 * exp(-b2 * x1)"
yvar <- paste("b1 * exp(-b2 * x1) * (1 + (1/", theta, ") * b1 * exp(-b2 * x1))", sep = "")
cfisher(ymean, yvar, ndpoints = 3, prec = 54)

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(LDOD)
Loading required package: Rsolnp
Loading required package: Rmpfr
Loading required package: gmp

Attaching package: 'gmp'

The following objects are masked from 'package:base':

    %*%, apply, crossprod, matrix, tcrossprod

C code of R package 'Rmpfr': GMP using 64 bits per limb


Attaching package: 'Rmpfr'

The following objects are masked from 'package:stats':

    dbinom, dnorm, dpois, pnorm

The following objects are masked from 'package:base':

    cbind, pmax, pmin, rbind

> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/LDOD/cfisher.Rd_%03d_medium.png", width=480, height=480)
> ### Name: cfisher
> ### Title: Auto-constructing Fisher Information matrix
> ### Aliases: cfisher
> ### Keywords: optimal design Fisher information matrix
> 
> ### ** Examples
> 
> ## Logistic dose response model
> ymean <- "(1/(exp(-b2 * (x1 - b1)) + 1))"
> yvar <- "(1/(exp(-b2 * (x1 - b1)) + 1)) * (1 - (1/(exp(-b2 * (x1 - b1)) + 1)))"
> res <- cfisher(ymean, yvar, ndpoints = 2, prec = 54)
> 
> # res$fim is Fisher information matrix for a two-points design 
> res$fim(x = c(x11 = 2, x21 = 3), w = c(w1 = .5, w2 = .5), b = c(b1 = .9, b2 = .8))
           [,1]       [,2]
[1,]  0.1086851 -0.2024087
[2,] -0.2024087  0.4173492
> 
> # res$fim is Fisher information matrix for a two-points design with 54 precision
> res$fim.mpfr(x = c(x11 = 2, x21 = 3), w = c(w1 = .5, w2 = .5), b = c(b1 = .9, b2 = .8))
'mpfrMatrix' of dim(.) =  (2, 2) of precision  54   bits 
     [,1]                  [,2]                 
[1,]  0.108685136448971514 -0.202408655633780798
[2,] -0.202408655633780798  0.417349208164616625
> 
> # Fisher information matrix for model:
> fim<- cfisher(ymean, yvar, ndpoints = 1, prec = 54)
> res$fim(x = c(x11 = 2), w = c(w1 = 1), b = c(b1 = .9, b2 = .8))
     [,1] [,2]
[1,]   NA   NA
[2,]   NA   NA
> 
> ## posison with E(y) = exp(b1 + b2 * x1 + b3 * x1^2 + b4 * x2 +b5 * x2^2 + b6 * x1 * x2)
> ymean <- yvar <- "exp(b1 + b2 * x1 + b3 * x1^2 + b4 * x2 +b5 * x2^2 + b6 * x1 * x2)"
> fim <- cfisher(ymean, yvar, ndpoints = 6, prec = 54)
> 
> # res$fim is Fisher information matrix for a six-points design 
> res$fim(x = c(1:12), w = rep(1/6, 6), b = c(1:6)) ## NAN
            [,1]        [,2]
[1,]  0.23666239 -0.03499786
[2,] -0.03499786  0.01749893
> 
> # res$fim.mpfr is Fisher information matrix for a six-points design with 53 precision
> res$fim.mpfr(x = c(1:12), w = rep(1/6, 6), b = c(1:6))
'mpfrMatrix' of dim(.) =  (2, 2) of precision  54   bits 
     [,1]                   [,2]                  
[1,]   0.236662390269004305 -0.0349978618011688274
[2,] -0.0349978618011688239  0.0174989309005844137
> 
> ## Linear regression with two indeoendent varibales (the design points are two-dimensional)
> ymean <- "b1 + b2 * x1 + b3 * x2"
> yvar = "1"
> res <- cfisher(ymean, yvar, ndpoints = 3, prec = 54)
> res$fim(x = c(1:6), w = c(.3, .3, .3))
     [,1] [,2] [,3]
[1,]  0.9  2.7  3.6
[2,]  2.7 10.5 13.2
[3,]  3.6 13.2 16.8
> res$fim.mpfr(x = c(1:6), w = c(.3, .3, .3))
'mpfrMatrix' of dim(.) =  (3, 3) of precision  54   bits 
     [,1]                 [,2]                [,3]               
[1,] 0.899999999999999967 2.69999999999999996 3.59999999999999987
[2,]  2.69999999999999996 10.5000000000000000 13.1999999999999993
[3,]  3.59999999999999987 13.1999999999999993 16.8000000000000007
> 
> ## Logistic model:
> ymean <- "1/(exp(-b1 - b2 * x1) + 1)"
> yvar <- "(1/(exp(-b1 - b2 * x1) + 1)) * (1 - (1/(exp(-b1 - b2 * x1) + 1)))"
> cfisher(ymean, yvar, ndpoints = 2, prec = 54)
$fim
function (x = c(x11, x21), w = c(w1, w2), b = c(b1, b2)) 
{
    Mat = matrix(c(+w[1] * 1/((1/(exp(-b[1] - b[2] * x[1]) + 
        1)) * (1 - (1/(exp(-b[1] - b[2] * x[1]) + 1)))) * exp(-b[1] - 
        b[2] * x[1])/(exp(-b[1] - b[2] * x[1]) + 1)^2 * exp(-b[1] - 
        b[2] * x[1])/(exp(-b[1] - b[2] * x[1]) + 1)^2 + w[2] * 
        1/((1/(exp(-b[1] - b[2] * x[2]) + 1)) * (1 - (1/(exp(-b[1] - 
        b[2] * x[2]) + 1)))) * exp(-b[1] - b[2] * x[2])/(exp(-b[1] - 
        b[2] * x[2]) + 1)^2 * exp(-b[1] - b[2] * x[2])/(exp(-b[1] - 
        b[2] * x[2]) + 1)^2, +w[1] * 1/((1/(exp(-b[1] - b[2] * 
        x[1]) + 1)) * (1 - (1/(exp(-b[1] - b[2] * x[1]) + 1)))) * 
        exp(-b[1] - b[2] * x[1]) * x[1]/(exp(-b[1] - b[2] * x[1]) + 
        1)^2 * exp(-b[1] - b[2] * x[1])/(exp(-b[1] - b[2] * x[1]) + 
        1)^2 + w[2] * 1/((1/(exp(-b[1] - b[2] * x[2]) + 1)) * 
        (1 - (1/(exp(-b[1] - b[2] * x[2]) + 1)))) * exp(-b[1] - 
        b[2] * x[2]) * x[2]/(exp(-b[1] - b[2] * x[2]) + 1)^2 * 
        exp(-b[1] - b[2] * x[2])/(exp(-b[1] - b[2] * x[2]) + 
        1)^2, +w[1] * 1/((1/(exp(-b[1] - b[2] * x[1]) + 1)) * 
        (1 - (1/(exp(-b[1] - b[2] * x[1]) + 1)))) * exp(-b[1] - 
        b[2] * x[1])/(exp(-b[1] - b[2] * x[1]) + 1)^2 * exp(-b[1] - 
        b[2] * x[1]) * x[1]/(exp(-b[1] - b[2] * x[1]) + 1)^2 + 
        w[2] * 1/((1/(exp(-b[1] - b[2] * x[2]) + 1)) * (1 - (1/(exp(-b[1] - 
            b[2] * x[2]) + 1)))) * exp(-b[1] - b[2] * x[2])/(exp(-b[1] - 
            b[2] * x[2]) + 1)^2 * exp(-b[1] - b[2] * x[2]) * 
            x[2]/(exp(-b[1] - b[2] * x[2]) + 1)^2, +w[1] * 1/((1/(exp(-b[1] - 
        b[2] * x[1]) + 1)) * (1 - (1/(exp(-b[1] - b[2] * x[1]) + 
        1)))) * exp(-b[1] - b[2] * x[1]) * x[1]/(exp(-b[1] - 
        b[2] * x[1]) + 1)^2 * exp(-b[1] - b[2] * x[1]) * x[1]/(exp(-b[1] - 
        b[2] * x[1]) + 1)^2 + w[2] * 1/((1/(exp(-b[1] - b[2] * 
        x[2]) + 1)) * (1 - (1/(exp(-b[1] - b[2] * x[2]) + 1)))) * 
        exp(-b[1] - b[2] * x[2]) * x[2]/(exp(-b[1] - b[2] * x[2]) + 
        1)^2 * exp(-b[1] - b[2] * x[2]) * x[2]/(exp(-b[1] - b[2] * 
        x[2]) + 1)^2), 2, 2)
    d = Mat
    return(d)
}
<environment: 0x4de0078>

$fim.mpfr
function (x = c(x11, x21), w = c(w1, w2), b = c(b1, b2)) 
{
    Mat = new("mpfrMatrix", c(+mpfr(w[1], precBits = 54) * 1/((1/(exp(-mpfr(b[1], 
        precBits = 54) - mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54)) + 1)) * (1 - (1/(exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[1], precBits = 54)) + 
        1)))) * exp(-mpfr(b[1], precBits = 54) - mpfr(b[2], precBits = 54) * 
        mpfr(x[1], precBits = 54))/(exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[1], precBits = 54)) + 
        1)^2 * exp(-mpfr(b[1], precBits = 54) - mpfr(b[2], precBits = 54) * 
        mpfr(x[1], precBits = 54))/(exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[1], precBits = 54)) + 
        1)^2 + mpfr(w[2], precBits = 54) * 1/((1/(exp(-mpfr(b[1], 
        precBits = 54) - mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54)) + 1)) * (1 - (1/(exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[2], precBits = 54)) + 
        1)))) * exp(-mpfr(b[1], precBits = 54) - mpfr(b[2], precBits = 54) * 
        mpfr(x[2], precBits = 54))/(exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[2], precBits = 54)) + 
        1)^2 * exp(-mpfr(b[1], precBits = 54) - mpfr(b[2], precBits = 54) * 
        mpfr(x[2], precBits = 54))/(exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[2], precBits = 54)) + 
        1)^2, +mpfr(w[1], precBits = 54) * 1/((1/(exp(-mpfr(b[1], 
        precBits = 54) - mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54)) + 1)) * (1 - (1/(exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[1], precBits = 54)) + 
        1)))) * exp(-mpfr(b[1], precBits = 54) - mpfr(b[2], precBits = 54) * 
        mpfr(x[1], precBits = 54)) * mpfr(x[1], precBits = 54)/(exp(-mpfr(b[1], 
        precBits = 54) - mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54)) + 1)^2 * exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[1], precBits = 54))/(exp(-mpfr(b[1], 
        precBits = 54) - mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54)) + 1)^2 + mpfr(w[2], precBits = 54) * 
        1/((1/(exp(-mpfr(b[1], precBits = 54) - mpfr(b[2], precBits = 54) * 
        mpfr(x[2], precBits = 54)) + 1)) * (1 - (1/(exp(-mpfr(b[1], 
        precBits = 54) - mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54)) + 1)))) * exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[2], precBits = 54)) * 
        mpfr(x[2], precBits = 54)/(exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[2], precBits = 54)) + 
        1)^2 * exp(-mpfr(b[1], precBits = 54) - mpfr(b[2], precBits = 54) * 
        mpfr(x[2], precBits = 54))/(exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[2], precBits = 54)) + 
        1)^2, +mpfr(w[1], precBits = 54) * 1/((1/(exp(-mpfr(b[1], 
        precBits = 54) - mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54)) + 1)) * (1 - (1/(exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[1], precBits = 54)) + 
        1)))) * exp(-mpfr(b[1], precBits = 54) - mpfr(b[2], precBits = 54) * 
        mpfr(x[1], precBits = 54))/(exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[1], precBits = 54)) + 
        1)^2 * exp(-mpfr(b[1], precBits = 54) - mpfr(b[2], precBits = 54) * 
        mpfr(x[1], precBits = 54)) * mpfr(x[1], precBits = 54)/(exp(-mpfr(b[1], 
        precBits = 54) - mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54)) + 1)^2 + mpfr(w[2], precBits = 54) * 
        1/((1/(exp(-mpfr(b[1], precBits = 54) - mpfr(b[2], precBits = 54) * 
        mpfr(x[2], precBits = 54)) + 1)) * (1 - (1/(exp(-mpfr(b[1], 
        precBits = 54) - mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54)) + 1)))) * exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[2], precBits = 54))/(exp(-mpfr(b[1], 
        precBits = 54) - mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54)) + 1)^2 * exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[2], precBits = 54)) * 
        mpfr(x[2], precBits = 54)/(exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[2], precBits = 54)) + 
        1)^2, +mpfr(w[1], precBits = 54) * 1/((1/(exp(-mpfr(b[1], 
        precBits = 54) - mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54)) + 1)) * (1 - (1/(exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[1], precBits = 54)) + 
        1)))) * exp(-mpfr(b[1], precBits = 54) - mpfr(b[2], precBits = 54) * 
        mpfr(x[1], precBits = 54)) * mpfr(x[1], precBits = 54)/(exp(-mpfr(b[1], 
        precBits = 54) - mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54)) + 1)^2 * exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[1], precBits = 54)) * 
        mpfr(x[1], precBits = 54)/(exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[1], precBits = 54)) + 
        1)^2 + mpfr(w[2], precBits = 54) * 1/((1/(exp(-mpfr(b[1], 
        precBits = 54) - mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54)) + 1)) * (1 - (1/(exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[2], precBits = 54)) + 
        1)))) * exp(-mpfr(b[1], precBits = 54) - mpfr(b[2], precBits = 54) * 
        mpfr(x[2], precBits = 54)) * mpfr(x[2], precBits = 54)/(exp(-mpfr(b[1], 
        precBits = 54) - mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54)) + 1)^2 * exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[2], precBits = 54)) * 
        mpfr(x[2], precBits = 54)/(exp(-mpfr(b[1], precBits = 54) - 
        mpfr(b[2], precBits = 54) * mpfr(x[2], precBits = 54)) + 
        1)^2), Dim = c(2L, 2L))
    d = Mat
    return(d)
}
<environment: 0x4de0078>

> 
> ## Poisson model:
> ymean <- yvar <-  "exp(b1 + b2 * x1)"
> cfisher(ymean, yvar, ndpoints = 2, prec = 54)
$fim
function (x = c(x11, x21), w = c(w1, w2), b = c(b1, b2)) 
{
    Mat = matrix(c(+w[1] * 1/(exp(b[1] + b[2] * x[1])) * exp(b[1] + 
        b[2] * x[1]) * exp(b[1] + b[2] * x[1]) + w[2] * 1/(exp(b[1] + 
        b[2] * x[2])) * exp(b[1] + b[2] * x[2]) * exp(b[1] + 
        b[2] * x[2]), +w[1] * 1/(exp(b[1] + b[2] * x[1])) * exp(b[1] + 
        b[2] * x[1]) * x[1] * exp(b[1] + b[2] * x[1]) + w[2] * 
        1/(exp(b[1] + b[2] * x[2])) * exp(b[1] + b[2] * x[2]) * 
        x[2] * exp(b[1] + b[2] * x[2]), +w[1] * 1/(exp(b[1] + 
        b[2] * x[1])) * exp(b[1] + b[2] * x[1]) * exp(b[1] + 
        b[2] * x[1]) * x[1] + w[2] * 1/(exp(b[1] + b[2] * x[2])) * 
        exp(b[1] + b[2] * x[2]) * exp(b[1] + b[2] * x[2]) * x[2], 
        +w[1] * 1/(exp(b[1] + b[2] * x[1])) * exp(b[1] + b[2] * 
            x[1]) * x[1] * exp(b[1] + b[2] * x[1]) * x[1] + w[2] * 
            1/(exp(b[1] + b[2] * x[2])) * exp(b[1] + b[2] * x[2]) * 
            x[2] * exp(b[1] + b[2] * x[2]) * x[2]), 2, 2)
    d = Mat
    return(d)
}
<environment: 0x4e8f5c0>

$fim.mpfr
function (x = c(x11, x21), w = c(w1, w2), b = c(b1, b2)) 
{
    Mat = new("mpfrMatrix", c(+mpfr(w[1], precBits = 54) * 1/(exp(mpfr(b[1], 
        precBits = 54) + mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54))) * exp(mpfr(b[1], precBits = 54) + mpfr(b[2], 
        precBits = 54) * mpfr(x[1], precBits = 54)) * exp(mpfr(b[1], 
        precBits = 54) + mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54)) + mpfr(w[2], precBits = 54) * 1/(exp(mpfr(b[1], 
        precBits = 54) + mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54))) * exp(mpfr(b[1], precBits = 54) + mpfr(b[2], 
        precBits = 54) * mpfr(x[2], precBits = 54)) * exp(mpfr(b[1], 
        precBits = 54) + mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54)), +mpfr(w[1], precBits = 54) * 1/(exp(mpfr(b[1], 
        precBits = 54) + mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54))) * exp(mpfr(b[1], precBits = 54) + mpfr(b[2], 
        precBits = 54) * mpfr(x[1], precBits = 54)) * mpfr(x[1], 
        precBits = 54) * exp(mpfr(b[1], precBits = 54) + mpfr(b[2], 
        precBits = 54) * mpfr(x[1], precBits = 54)) + mpfr(w[2], 
        precBits = 54) * 1/(exp(mpfr(b[1], precBits = 54) + mpfr(b[2], 
        precBits = 54) * mpfr(x[2], precBits = 54))) * exp(mpfr(b[1], 
        precBits = 54) + mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54)) * mpfr(x[2], precBits = 54) * exp(mpfr(b[1], 
        precBits = 54) + mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54)), +mpfr(w[1], precBits = 54) * 1/(exp(mpfr(b[1], 
        precBits = 54) + mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54))) * exp(mpfr(b[1], precBits = 54) + mpfr(b[2], 
        precBits = 54) * mpfr(x[1], precBits = 54)) * exp(mpfr(b[1], 
        precBits = 54) + mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54)) * mpfr(x[1], precBits = 54) + mpfr(w[2], 
        precBits = 54) * 1/(exp(mpfr(b[1], precBits = 54) + mpfr(b[2], 
        precBits = 54) * mpfr(x[2], precBits = 54))) * exp(mpfr(b[1], 
        precBits = 54) + mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54)) * exp(mpfr(b[1], precBits = 54) + mpfr(b[2], 
        precBits = 54) * mpfr(x[2], precBits = 54)) * mpfr(x[2], 
        precBits = 54), +mpfr(w[1], precBits = 54) * 1/(exp(mpfr(b[1], 
        precBits = 54) + mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54))) * exp(mpfr(b[1], precBits = 54) + mpfr(b[2], 
        precBits = 54) * mpfr(x[1], precBits = 54)) * mpfr(x[1], 
        precBits = 54) * exp(mpfr(b[1], precBits = 54) + mpfr(b[2], 
        precBits = 54) * mpfr(x[1], precBits = 54)) * mpfr(x[1], 
        precBits = 54) + mpfr(w[2], precBits = 54) * 1/(exp(mpfr(b[1], 
        precBits = 54) + mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54))) * exp(mpfr(b[1], precBits = 54) + mpfr(b[2], 
        precBits = 54) * mpfr(x[2], precBits = 54)) * mpfr(x[2], 
        precBits = 54) * exp(mpfr(b[1], precBits = 54) + mpfr(b[2], 
        precBits = 54) * mpfr(x[2], precBits = 54)) * mpfr(x[2], 
        precBits = 54)), Dim = c(2L, 2L))
    d = Mat
    return(d)
}
<environment: 0x4e8f5c0>

> 
> ## Poisson dose response model:
> ymean <- yvar <- "b1 * exp(-b2 * x1)"
> cfisher(ymean, yvar, ndpoints = 2, prec = 54)
$fim
function (x = c(x11, x21), w = c(w1, w2), b = c(b1, b2)) 
{
    Mat = matrix(c(+w[1] * 1/(b[1] * exp(-b[2] * x[1])) * exp(-b[2] * 
        x[1]) * exp(-b[2] * x[1]) + w[2] * 1/(b[1] * exp(-b[2] * 
        x[2])) * exp(-b[2] * x[2]) * exp(-b[2] * x[2]), +w[1] * 
        1/(b[1] * exp(-b[2] * x[1])) * -(b[1] * (exp(-b[2] * 
        x[1]) * x[1])) * exp(-b[2] * x[1]) + w[2] * 1/(b[1] * 
        exp(-b[2] * x[2])) * -(b[1] * (exp(-b[2] * x[2]) * x[2])) * 
        exp(-b[2] * x[2]), +w[1] * 1/(b[1] * exp(-b[2] * x[1])) * 
        exp(-b[2] * x[1]) * -(b[1] * (exp(-b[2] * x[1]) * x[1])) + 
        w[2] * 1/(b[1] * exp(-b[2] * x[2])) * exp(-b[2] * x[2]) * 
            -(b[1] * (exp(-b[2] * x[2]) * x[2])), +w[1] * 1/(b[1] * 
        exp(-b[2] * x[1])) * -(b[1] * (exp(-b[2] * x[1]) * x[1])) * 
        -(b[1] * (exp(-b[2] * x[1]) * x[1])) + w[2] * 1/(b[1] * 
        exp(-b[2] * x[2])) * -(b[1] * (exp(-b[2] * x[2]) * x[2])) * 
        -(b[1] * (exp(-b[2] * x[2]) * x[2]))), 2, 2)
    d = Mat
    return(d)
}
<environment: 0x4f0ed50>

$fim.mpfr
function (x = c(x11, x21), w = c(w1, w2), b = c(b1, b2)) 
{
    Mat = new("mpfrMatrix", c(+mpfr(w[1], precBits = 54) * 1/(mpfr(b[1], 
        precBits = 54) * exp(-mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54))) * exp(-mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54)) * exp(-mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54)) + mpfr(w[2], precBits = 54) * 1/(mpfr(b[1], 
        precBits = 54) * exp(-mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54))) * exp(-mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54)) * exp(-mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54)), +mpfr(w[1], precBits = 54) * 1/(mpfr(b[1], 
        precBits = 54) * exp(-mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54))) * -(mpfr(b[1], precBits = 54) * (exp(-mpfr(b[2], 
        precBits = 54) * mpfr(x[1], precBits = 54)) * mpfr(x[1], 
        precBits = 54))) * exp(-mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54)) + mpfr(w[2], precBits = 54) * 1/(mpfr(b[1], 
        precBits = 54) * exp(-mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54))) * -(mpfr(b[1], precBits = 54) * (exp(-mpfr(b[2], 
        precBits = 54) * mpfr(x[2], precBits = 54)) * mpfr(x[2], 
        precBits = 54))) * exp(-mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54)), +mpfr(w[1], precBits = 54) * 1/(mpfr(b[1], 
        precBits = 54) * exp(-mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54))) * exp(-mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54)) * -(mpfr(b[1], precBits = 54) * (exp(-mpfr(b[2], 
        precBits = 54) * mpfr(x[1], precBits = 54)) * mpfr(x[1], 
        precBits = 54))) + mpfr(w[2], precBits = 54) * 1/(mpfr(b[1], 
        precBits = 54) * exp(-mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54))) * exp(-mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54)) * -(mpfr(b[1], precBits = 54) * (exp(-mpfr(b[2], 
        precBits = 54) * mpfr(x[2], precBits = 54)) * mpfr(x[2], 
        precBits = 54))), +mpfr(w[1], precBits = 54) * 1/(mpfr(b[1], 
        precBits = 54) * exp(-mpfr(b[2], precBits = 54) * mpfr(x[1], 
        precBits = 54))) * -(mpfr(b[1], precBits = 54) * (exp(-mpfr(b[2], 
        precBits = 54) * mpfr(x[1], precBits = 54)) * mpfr(x[1], 
        precBits = 54))) * -(mpfr(b[1], precBits = 54) * (exp(-mpfr(b[2], 
        precBits = 54) * mpfr(x[1], precBits = 54)) * mpfr(x[1], 
        precBits = 54))) + mpfr(w[2], precBits = 54) * 1/(mpfr(b[1], 
        precBits = 54) * exp(-mpfr(b[2], precBits = 54) * mpfr(x[2], 
        precBits = 54))) * -(mpfr(b[1], precBits = 54) * (exp(-mpfr(b[2], 
        precBits = 54) * mpfr(x[2], precBits = 54)) * mpfr(x[2], 
        precBits = 54))) * -(mpfr(b[1], precBits = 54) * (exp(-mpfr(b[2], 
        precBits = 54) * mpfr(x[2], precBits = 54)) * mpfr(x[2], 
        precBits = 54)))), Dim = c(2L, 2L))
    d = Mat
    return(d)
}
<environment: 0x4f0ed50>

> 
> ## Inverse Quadratic model:
> ymean <- "(b1 * x1)/(b2 + x1 + b3 * (x1)^2)"
> yvar <- "1"
> cfisher(ymean, yvar, ndpoints = 3, prec = 54)
$fim
function (x = c(x11, x21, x31), w = c(w1, w2, w3), b = c(b1, 
    b2, b3)) 
{
    Mat = matrix(c(+w[1] * 1/(1) * x[1]/(b[2] + x[1] + b[3] * 
        (x[1])^2) * x[1]/(b[2] + x[1] + b[3] * (x[1])^2) + w[2] * 
        1/(1) * x[2]/(b[2] + x[2] + b[3] * (x[2])^2) * x[2]/(b[2] + 
        x[2] + b[3] * (x[2])^2) + w[3] * 1/(1) * x[3]/(b[2] + 
        x[3] + b[3] * (x[3])^2) * x[3]/(b[2] + x[3] + b[3] * 
        (x[3])^2), +w[1] * 1/(1) * -((b[1] * x[1])/(b[2] + x[1] + 
        b[3] * (x[1])^2)^2) * x[1]/(b[2] + x[1] + b[3] * (x[1])^2) + 
        w[2] * 1/(1) * -((b[1] * x[2])/(b[2] + x[2] + b[3] * 
            (x[2])^2)^2) * x[2]/(b[2] + x[2] + b[3] * (x[2])^2) + 
        w[3] * 1/(1) * -((b[1] * x[3])/(b[2] + x[3] + b[3] * 
            (x[3])^2)^2) * x[3]/(b[2] + x[3] + b[3] * (x[3])^2), 
        +w[1] * 1/(1) * -((b[1] * x[1]) * (x[1])^2/(b[2] + x[1] + 
            b[3] * (x[1])^2)^2) * x[1]/(b[2] + x[1] + b[3] * 
            (x[1])^2) + w[2] * 1/(1) * -((b[1] * x[2]) * (x[2])^2/(b[2] + 
            x[2] + b[3] * (x[2])^2)^2) * x[2]/(b[2] + x[2] + 
            b[3] * (x[2])^2) + w[3] * 1/(1) * -((b[1] * x[3]) * 
            (x[3])^2/(b[2] + x[3] + b[3] * (x[3])^2)^2) * x[3]/(b[2] + 
            x[3] + b[3] * (x[3])^2), +w[1] * 1/(1) * x[1]/(b[2] + 
            x[1] + b[3] * (x[1])^2) * -((b[1] * x[1])/(b[2] + 
            x[1] + b[3] * (x[1])^2)^2) + w[2] * 1/(1) * x[2]/(b[2] + 
            x[2] + b[3] * (x[2])^2) * -((b[1] * x[2])/(b[2] + 
            x[2] + b[3] * (x[2])^2)^2) + w[3] * 1/(1) * x[3]/(b[2] + 
            x[3] + b[3] * (x[3])^2) * -((b[1] * x[3])/(b[2] + 
            x[3] + b[3] * (x[3])^2)^2), +w[1] * 1/(1) * -((b[1] * 
            x[1])/(b[2] + x[1] + b[3] * (x[1])^2)^2) * -((b[1] * 
            x[1])/(b[2] + x[1] + b[3] * (x[1])^2)^2) + w[2] * 
            1/(1) * -((b[1] * x[2])/(b[2] + x[2] + b[3] * (x[2])^2)^2) * 
            -((b[1] * x[2])/(b[2] + x[2] + b[3] * (x[2])^2)^2) + 
            w[3] * 1/(1) * -((b[1] * x[3])/(b[2] + x[3] + b[3] * 
                (x[3])^2)^2) * -((b[1] * x[3])/(b[2] + x[3] + 
                b[3] * (x[3])^2)^2), +w[1] * 1/(1) * -((b[1] * 
            x[1]) * (x[1])^2/(b[2] + x[1] + b[3] * (x[1])^2)^2) * 
            -((b[1] * x[1])/(b[2] + x[1] + b[3] * (x[1])^2)^2) + 
            w[2] * 1/(1) * -((b[1] * x[2]) * (x[2])^2/(b[2] + 
                x[2] + b[3] * (x[2])^2)^2) * -((b[1] * x[2])/(b[2] + 
                x[2] + b[3] * (x[2])^2)^2) + w[3] * 1/(1) * -((b[1] * 
            x[3]) * (x[3])^2/(b[2] + x[3] + b[3] * (x[3])^2)^2) * 
            -((b[1] * x[3])/(b[2] + x[3] + b[3] * (x[3])^2)^2), 
        +w[1] * 1/(1) * x[1]/(b[2] + x[1] + b[3] * (x[1])^2) * 
            -((b[1] * x[1]) * (x[1])^2/(b[2] + x[1] + b[3] * 
                (x[1])^2)^2) + w[2] * 1/(1) * x[2]/(b[2] + x[2] + 
            b[3] * (x[2])^2) * -((b[1] * x[2]) * (x[2])^2/(b[2] + 
            x[2] + b[3] * (x[2])^2)^2) + w[3] * 1/(1) * x[3]/(b[2] + 
            x[3] + b[3] * (x[3])^2) * -((b[1] * x[3]) * (x[3])^2/(b[2] + 
            x[3] + b[3] * (x[3])^2)^2), +w[1] * 1/(1) * -((b[1] * 
            x[1])/(b[2] + x[1] + b[3] * (x[1])^2)^2) * -((b[1] * 
            x[1]) * (x[1])^2/(b[2] + x[1] + b[3] * (x[1])^2)^2) + 
            w[2] * 1/(1) * -((b[1] * x[2])/(b[2] + x[2] + b[3] * 
                (x[2])^2)^2) * -((b[1] * x[2]) * (x[2])^2/(b[2] + 
                x[2] + b[3] * (x[2])^2)^2) + w[3] * 1/(1) * -((b[1] * 
            x[3])/(b[2] + x[3] + b[3] * (x[3])^2)^2) * -((b[1] * 
            x[3]) * (x[3])^2/(b[2] + x[3] + b[3] * (x[3])^2)^2), 
        +w[1] * 1/(1) * -((b[1] * x[1]) * (x[1])^2/(b[2] + x[1] + 
            b[3] * (x[1])^2)^2) * -((b[1] * x[1]) * (x[1])^2/(b[2] + 
            x[1] + b[3] * (x[1])^2)^2) + w[2] * 1/(1) * -((b[1] * 
            x[2]) * (x[2])^2/(b[2] + x[2] + b[3] * (x[2])^2)^2) * 
            -((b[1] * x[2]) * (x[2])^2/(b[2] + x[2] + b[3] * 
                (x[2])^2)^2) + w[3] * 1/(1) * -((b[1] * x[3]) * 
            (x[3])^2/(b[2] + x[3] + b[3] * (x[3])^2)^2) * -((b[1] * 
            x[3]) * (x[3])^2/(b[2] + x[3] + b[3] * (x[3])^2)^2)), 
        3, 3)
    d = Mat
    return(d)
}
<environment: 0x4f8f9f0>

$fim.mpfr
function (x = c(x11, x21, x31), w = c(w1, w2, w3), b = c(b1, 
    b2, b3)) 
{
    Mat = new("mpfrMatrix", c(+mpfr(w[1], precBits = 54) * 1/(1) * 
        mpfr(x[1], precBits = 54)/(mpfr(b[2], precBits = 54) + 
        mpfr(x[1], precBits = 54) + mpfr(b[3], precBits = 54) * 
        (mpfr(x[1], precBits = 54))^2) * mpfr(x[1], precBits = 54)/(mpfr(b[2], 
        precBits = 54) + mpfr(x[1], precBits = 54) + mpfr(b[3], 
        precBits = 54) * (mpfr(x[1], precBits = 54))^2) + mpfr(w[2], 
        precBits = 54) * 1/(1) * mpfr(x[2], precBits = 54)/(mpfr(b[2], 
        precBits = 54) + mpfr(x[2], precBits = 54) + mpfr(b[3], 
        precBits = 54) * (mpfr(x[2], precBits = 54))^2) * mpfr(x[2], 
        precBits = 54)/(mpfr(b[2], precBits = 54) + mpfr(x[2], 
        precBits = 54) + mpfr(b[3], precBits = 54) * (mpfr(x[2], 
        precBits = 54))^2) + mpfr(w[3], precBits = 54) * 1/(1) * 
        mpfr(x[3], precBits = 54)/(mpfr(b[2], precBits = 54) + 
        mpfr(x[3], precBits = 54) + mpfr(b[3], precBits = 54) * 
        (mpfr(x[3], precBits = 54))^2) * mpfr(x[3], precBits = 54)/(mpfr(b[2], 
        precBits = 54) + mpfr(x[3], precBits = 54) + mpfr(b[3], 
        precBits = 54) * (mpfr(x[3], precBits = 54))^2), +mpfr(w[1], 
        precBits = 54) * 1/(1) * -((mpfr(b[1], precBits = 54) * 
        mpfr(x[1], precBits = 54))/(mpfr(b[2], precBits = 54) + 
        mpfr(x[1], precBits = 54) + mpfr(b[3], precBits = 54) * 
        (mpfr(x[1], precBits = 54))^2)^2) * mpfr(x[1], precBits = 54)/(mpfr(b[2], 
        precBits = 54) + mpfr(x[1], precBits = 54) + mpfr(b[3], 
        precBits = 54) * (mpfr(x[1], precBits = 54))^2) + mpfr(w[2], 
        precBits = 54) * 1/(1) * -((mpfr(b[1], precBits = 54) * 
        mpfr(x[2], precBits = 54))/(mpfr(b[2], precBits = 54) + 
        mpfr(x[2], precBits = 54) + mpfr(b[3], precBits = 54) * 
        (mpfr(x[2], precBits = 54))^2)^2) * mpfr(x[2], precBits = 54)/(mpfr(b[2], 
        precBits = 54) + mpfr(x[2], precBits = 54) + mpfr(b[3], 
        precBits = 54) * (mpfr(x[2], precBits = 54))^2) + mpfr(w[3], 
        precBits = 54) * 1/(1) * -((mpfr(b[1], precBits = 54) * 
        mpfr(x[3], precBits = 54))/(mpfr(b[2], precBits = 54) + 
        mpfr(x[3], precBits = 54) + mpfr(b[3], precBits = 54) * 
        (mpfr(x[3], precBits = 54))^2)^2) * mpfr(x[3], precBits = 54)/(mpfr(b[2], 
        precBits = 54) + mpfr(x[3], precBits = 54) + mpfr(b[3], 
        precBits = 54) * (mpfr(x[3], precBits = 54))^2), +mpfr(w[1], 
        precBits = 54) * 1/(1) * -((mpfr(b[1], precBits = 54) * 
        mpfr(x[1], precBits = 54)) * (mpfr(x[1], precBits = 54))^2/(mpfr(b[2], 
        precBits = 54) + mpfr(x[1], precBits = 54) + mpfr(b[3], 
        precBits = 54) * (mpfr(x[1], precBits = 54))^2)^2) * 
        mpfr(x[1], precBits = 54)/(mpfr(b[2], precBits = 54) + 
        mpfr(x[1], precBits = 54) + mpfr(b[3], precBits = 54) * 
        (mpfr(x[1], precBits = 54))^2) + mpfr(w[2], precBits = 54) * 
        1/(1) * -((mpfr(b[1], precBits = 54) * mpfr(x[2], precBits = 54)) * 
        (mpfr(x[2], precBits = 54))^2/(mpfr(b[2], precBits = 54) + 
        mpfr(x[2], precBits = 54) + mpfr(b[3], precBits = 54) * 
        (mpfr(x[2], precBits = 54))^2)^2) * mpfr(x[2], precBits = 54)/(mpfr(b[2], 
        precBits = 54) + mpfr(x[2], precBits = 54) + mpfr(b[3], 
        precBits = 54) * (mpfr(x[2], precBits = 54))^2) + mpfr(w[3], 
        precBits = 54) * 1/(1) * -((mpfr(b[1], precBits = 54) * 
        mpfr(x[3], precBits = 54)) * (mpfr(x[3], precBits = 54))^2/(mpfr(b[2], 
        precBits = 54) + mpfr(x[3], precBits = 54) + mpfr(b[3], 
        precBits = 54) * (mpfr(x[3], precBits = 54))^2)^2) * 
        mpfr(x[3], precBits = 54)/(mpfr(b[2], precBits = 54) + 
        mpfr(x[3], precBits = 54) + mpfr(b[3], precBits = 54) * 
        (mpfr(x[3], precBits = 54))^2), +mpfr(w[1], precBits = 54) * 
        1/(1) * mpfr(x[1], precBits = 54)/(mpfr(b[2], precBits = 54) + 
        mpfr(x[1], precBits = 54) + mpfr(b[3], precBits = 54) * 
        (mpfr(x[1], precBits = 54))^2) * -((mpfr(b[1], precBits = 54) * 
        mpfr(x[1], precBits = 54))/(mpfr(b[2], precBits = 54) + 
        mpfr(x[1], precBits = 54) + mpfr(b[3], precBits = 54) * 
        (mpfr(x[1], precBits = 54))^2)^2) + mpfr(w[2], precBits = 54) * 
        1/(1) * mpfr(x[2], precBits = 54)/(mpfr(b[2], precBits = 54) + 
        mpfr(x[2], precBits = 54) + mpfr(b[3], precBits = 54) * 
        (mpfr(x[2], precBits = 54))^2) * -((mpfr(b[1], precBits = 54) * 
        mpfr(x[2], precBits = 54))/(mpfr(b[2], precBits = 54) + 
        mpfr(x[2], precBits = 54) + mpfr(b[3], precBits = 54) * 
        (mpfr(x[2], precBits = 54))^2)^2) + mpfr(w[3], precBits = 54) * 
        1/(1) * mpfr(x[3], precBits = 54)/(mpfr(b[2], precBits = 54) + 
        mpfr(x[3], precBits = 54) + mpfr(b[3], precBits = 54) * 
        (mpfr(x[3], precBits = 54))^2) * -((mpfr(b[1