text:
# using importance sampling to reduce the error of outscattered radiance

http://renderman.pixar.com/view/implementing-a-skin-bssrdf

In this paper on subsurface scattering, I'm trying to understand how importance sampling is used to compute single scattering.

It says that the outscattered radiance for the single term, $L_o$, can be rewritten as the product of the exponential falloff in $s'_o$ with another function:

$L_o(x_o,\vec{\omega}_o) = \bigg[ L_i(x_i,\vec{\omega}_i)\cfrac{\sigma_s(x_o)F_p(\vec{\omega}'_i \cdotp \vec{\omega}'_o )}{\sigma_{tc}}e^{-s'_i \sigma_t(x_i)} \bigg] e^{-s_o\sigma(x_o)}$

By picking samples according to $x \sim \sigma_t e^{-\sigma_t x}$, the integration in $L_o$ can be approximated by a summation by applying importance sampling. Also by choosing $s'_o = \cfrac{-\mathbb{log(random())}}{\sigma_{tr}}$, you can sum up all contributions without computing the falloff term in $s_o$.

I understand you can approximate an integral using importance sampling, but I didn't follow how it applies to the rendering equation. Please could you describe how this works in greater detail?

$E[f(x)]=\int f(x)p(x) = \int f(x)\cfrac{p(x)}{q(x)}q(x) \approx \cfrac{1}{n}\sum \limits_{i=1}^n f(x_i)\cfrac{p(x_i)}{q(x_i)}$

Also, where does single scattering fit into the general rendering equation?

$L_{\text{o}}(\mathbf x,\, \omega_{\text{o}}) \,=\, L_e(\mathbf x,\, \omega_{\text{o}}) \ +\, \int_\Omega f_r(\mathbf x,\, \omega_{\text{i}},\, \omega_{\text{o}})\, L_{\text{i}}(\mathbf x,\, \omega_{\text{i}})\, (\omega_{\text{i}}\,\cdot\,\mathbf n)\, \operatorname d \omega_{\text{i}}$

In general, the idea of importance sampling is to distribute samples in a way that matches the function being integrated (or more practically, matching some factor in it). This reduces the variance of the output values, allowing the Monte Carlo integration to converge faster.

In volume scattering, when a medium is homogeneous, we have the Beer–Lambert law, which states that the transmittance falls off exponentially with depth. That's the source of factors like $e^{-\sigma x}$ in those equations. So when raymarching through a medium, the regions closer to the start of the ray generally contribute more light to the output, while the farther regions have more of their light absorbed along the way.

So, if you sample points along the ray from an exponential distribution, the samples cluster toward the closer, more important regions. Then you just have to evaluate the scattering at those points and average over them.

BTW, importance sampling in the context of surface scattering is explained pretty well by this answer to another question.

On your second question, "where does single scattering fit into the general rendering equation?", well, it doesn't! What people call the "general rendering equation" actually only includes surface scattering, not volume scattering. So it's not really that "general". A general rendering equation for volume scattering might look something like:

$$\nabla_\omega \, L(x, \omega) = \sigma_e(x) L_e(x, \omega) - \sigma_t(x) L(x, \omega) + \sigma_s(x) \int p(\omega, \omega') \, L(x, \omega') \, d\omega'$$

This states that the derivative of radiance as you move along a ray is given by the sum of emitted light from the volume, minus extinction, plus in-scattering.

To combine surface and volume scattering, we usually just say that the surface version of the equation applies at points on surfaces, and the volume version applies everywhere else.

• Is the dipole approximation at all related to the Beer–Lambert law? Jan 22 '16 at 3:16
• @maogenc The dipole approximation concerns homogeneous materials with a very high scattering coefficient, such that they're dominated by multiple scattering. The Beer–Lambert law would also apply to such materials. Jan 22 '16 at 6:38


text:
## Regression & Finite Mixture Models

I wrote a post a while back about Mixture Distributions and Model Comparisons. This post continues on that theme and tries to model multiple data generating processes into a single model. The code for this post is available at the github repository. There were many useful resources that helped me understand this model, and some are mentioned below and others are mentioned in the previous post in this topic:

I use the R library flexmix, for the sample data and fitting the model, in addition to Stan and MCMC approach. The data distribution is shown in the figure below, where looking at the figure suggests that there may be two data generation processes at work here. To start with, we try and find a distribution that will capture certain aspects of this data, and use posterior predictive checks to see if our distribution choice is appropriate. A previous post on posterior predictive checks covers this topic in more detail, and tests different distributions. To keep it short, I just try a normal mixture distribution with 2 mixture components and then perform some tests.

The figure in the left panel shows the density plot for the raw data, and the panel on the right shows the same data, with smattering of 200 predictive data sets from the fitted distribution – visually it appears to be fine, with the left component slightly higher in some cases. However the aspects of the data that I want to our distribution to capture are the mean, variance, maximum and minimum values, and the test result PValues suggest that any variation can be explained by sampling variability.

```> ivTestQuantities
variance  minimum  maximum     mean
0.380    0.155    0.300    0.420
```

The first steps are to fit a mixture regression model with 2 components and one predictor using the flexmix package and then to repeat this using Stan.

```> fit.flex = flexmix(yn ~ x, data=NPreg, k=2)
> summary(fit.flex)

Call:
flexmix(formula = yn ~ x, data = NPreg, k = 2)

prior size post>0 ratio
Comp.1 0.561  104    198 0.525
Comp.2 0.439   96    146 0.658

'log Lik.' -704.0635 (df=7)
AIC: 1422.127   BIC: 1445.215

> ## fitted coefficients
> parameters(fit.flex)
Comp.1     Comp.2
coef.(Intercept) 31.7644433 -0.7685309
coef.x            0.1867837  5.1841287
sigma             7.7718438  3.3828795
```

The same model is now broken down into parts and fit using Stan – you can write a log posterior function in R as well to do this, but I find that the optimiser tends to not explore the parameter space if we track the mixture components as well. However if we keep the mixture components fixed, then it tends to converge to the same point as Stan.

```> print(fit.stan, digi=3)
Inference for Stan model: fitNormalMixtureRegression.
1 chains, each with iter=4000; warmup=2000; thin=1;
post-warmup draws per chain=2000, total post-warmup draws=2000.

mean se_mean    sd     2.5%      25%      50%      75%    97.5% n_eff  Rhat
sigma[1]          3.405   0.010 0.395    2.676    3.126    3.395    3.667    4.223  1448 1.002
sigma[2]          7.848   0.019 0.632    6.725    7.407    7.805    8.270    9.164  1142 1.000
iMixWeights[1]    0.439   0.001 0.046    0.350    0.406    0.438    0.472    0.530  1000 1.000
iMixWeights[2]    0.561   0.001 0.046    0.470    0.528    0.562    0.594    0.650  1000 1.000
betasMix1[1]      5.163   0.003 0.146    4.870    5.067    5.160    5.262    5.452  2000 1.001
betasMix2[1]      0.176   0.013 0.347   -0.501   -0.056    0.171    0.391    0.861   751 1.001
mu[1]            -0.686   0.016 0.713   -2.072   -1.164   -0.705   -0.228    0.750  2000 1.000
mu[2]            31.846   0.067 1.971   27.898   30.590   31.890   33.136   35.791   853 1.000
lp__           -708.099   0.063 1.885 -712.367 -709.142 -707.715 -706.712 -705.449   899 1.001

Samples were drawn using NUTS(diag_e) at Tue Jul  4 16:10:02 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
```

Comparing the results from the two procedures, i.e. flexmix and the Stan model, suggest a nice convergence between the two methods.

## Fitted or Predicted Values

Once the coefficients of the regression equation have been determined, we can calculate the predicted (on new data) or fitted (current data) values. There are two ways to get this data: either get predicted data for each component; Or get an aggregate predicted value using each component.

```> intercepts
[1] -0.6598333 31.5694812
> betas
[1] 5.1614796 0.2241089
> iMixWeights
[1] 0.44 0.56
> iPred1 = lStanData\$X %*% c(intercepts[1], betas[1])
> iPred2 = lStanData\$X %*% c(intercepts[2], betas[2])
> ## compare with fit.flex
> head(f)
Comp.1    Comp.2
[1,] 32.54457 20.883670
[2,] 31.98889  5.460877
[3,] 32.19311 11.129077
[4,] 32.87877 30.159296
[5,] 32.20489 11.456077
[6,] 33.19092 38.822976
> head(cbind(iPred1, iPred2))
[,1]     [,2]
1 20.897771 32.50550
2  5.542358 31.83878
3 11.185795 32.08381
4 30.132872 32.90649
5 11.511366 32.09795
6 38.758702 33.28101
```

However getting one aggregate predicted value for the response variable may be more appropriate in certain settings, and will require us to take a weighted average of each component (the weight determined by the mixing probability).

```> iAggregate = cbind(iPred1, iPred2)
> iAggregate = sweep(iAggregate, 2, iMixWeights, '*')
> iAggregate = rowSums(iAggregate)
> dfPredicted = data.frame(stan=iAggregate, flexmix=p.agg)
> head(dfPredicted)
stan  flexmix
1 27.39810 27.42164
2 20.26835 20.33447
3 22.88868 22.93915
4 31.68610 31.68404
5 23.03985 23.08942
6 35.69120 35.66522
> ## calculate Mean Squared Error MSE
> mean((dfPredicted\$flexmix - NPreg\$yn)^2)
[1] 104.4622
> mean((dfPredicted\$stan - NPreg\$yn)^2)
[1] 104.3325
> ## both pretty close
```

Comparing the aggregate values from the two models show that both have produced similar results, and the mean squared errors in each case are also similar.

## Model Checking: Scoring and Comparing Models

This is another post in the series of model checking posts. Previously we looked at which aspects of the data and model are compatible, using posterior predictive checks. Once we have selected a model or a set of models for the data, we would like to score and compare them. One aspect of comparison using Mixture Distributions and Bayes Factors has been show in a previous post. Here we use posterior predictive checks to calculate information criteria and leave one out cross validation (and adjusting for number of parameters in the model). The examples and algorithms used in the post can be seen in more detail in:

[1] Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B. (2013). Bayesian Data Analysis, Third Edition (Texts in Statistical Science). Book.

[2] Albert, J., Gentleman, R., Parmigiani, G., & Hornik, K. (2009). Bayesian computation with R. Bayesian Computation with R. http://doi.org/10.1007/978-0-387-92298-0

[3] James, G., Witten, D., Hastie, T., & Tibshirani, R. (2013). An introduction to statistical learning: with applications in R. Springer texts in statistics (Vol. XIV). http://doi.org/10.1007/978-1-4614-7138-7

The predictive accuracy of a model can be used to evaluate and compare a model against other candidate models. The accuracy can be defined using some form of a function, e.g. mean squared error, log-likelihood, log predictive density etc. This ties into the concept of deviance or information criteria, which puts different models (with different parameters) on a common scale (along with some adjustments made for the number of parameters being estimated).

The data to calculate this predictive accuracy can be recycled (within sample) or external (out of sample i.e. new data) or using cross-validation (Ref 3 has a good chapter on cross-validation).

## Example Data and Model

The data used in this example is from the ‘bread and peace’ model, and the paper can be seen here, while the example is from Ref [1], and the R code can be found in the github repository.

We first fit a linear model using the lm function in R.

```Call:
lm(formula = VoteShare ~ IncomeGrowth, data = dfData)

Residuals:
Min      1Q  Median      3Q     Max
-8.8370 -0.3748  0.1379  1.7745  5.6291

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   45.8027     1.7297  26.480 1.07e-12 ***
IncomeGrowth   3.1809     0.7226   4.402 0.000715 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.816 on 13 degrees of freedom
Multiple R-squared:  0.5985,    Adjusted R-squared:  0.5676
F-statistic: 19.38 on 1 and 13 DF,  p-value: 0.0007149
```

There are three parameters being estimated here, the two coefficients and the standard deviation for the likelihood function.

```## write the model with the log posterior function
# use the likelihood with non-informative priors
mylogpost = function(theta, data){
sigma = exp(theta['sigma'])
betas = theta[-1]
x = data\$pred
y = data\$resp

# get the predicted or fitted value
mModMatrix = model.matrix(y ~ x)
mCoef = matrix(betas, nrow = length(betas))
iFitted = mModMatrix %*% mCoef # using identity link
## likelihood function
llik = sum(dnorm(y, iFitted, sigma, log=T))
## prior, non informative
lprior = 1
lpost = llik + lprior
return(lpost)
}
```

We can define this regression model using our custom function mylogpost and we are tracking 3 parameters (sigma, and 2 betas) – we are also using non-informative priors so our results should be similar to the lm model.

```> fit = laplace(mylogpost, start, lData)
> fit
\$mode
sigma     beta0     beta1
1.268109 45.806382  3.180065

\$var
sigma         beta0         beta1
sigma  3.335900e-02  0.0002461944 -5.857976e-05
beta0  2.461944e-04  2.5949267841 -8.909553e-01
beta1 -5.857976e-05 -0.8909553156  4.528745e-01

\$int
[1] -38.72532

\$converge
[1] TRUE

> se = sqrt(diag(fit\$var))
> se
sigma     beta0     beta1
0.1826445 1.6108776 0.6729595
## coef and sd
> fit\$mode[-1]
beta0     beta1
45.806382  3.180065
> exp(fit\$mode[1])
sigma
3.554124
>
```

We take a sample from the posterior distribution of these three parameters using a multivariate t proposal density, while the target density is mylogpost using a Sampling Importance Resampling (SIR) algorithm (I will write a post about these sampling methods at some point in the future). You can try using the approach suggested in the post in parameter estimation, but if these are sampled independently then we can not account for the covariance between the parameters – which can be avoided by using a multivariate normal sampling.

```## parameters for the multivariate t density
tpar = list(m=fit\$mode, var=fit\$var*2, df=4)
## get a sample directly and using sir (sampling importance resampling with a t proposal density)
s = sir(mylogpost, tpar, 1000, lData)
#s = rmt(10000, fit\$mode, S = fit\$var)
sigSample = s[,'sigma']
beta0Sample = s[,'beta0']
beta1Sample = s[,'beta1']
```

We use the SIR algorithm from [2], but you can write your own, its straight forward once you know how. The tpar is a list of parameters for the multivariate t proposal density. Just to show how our sampler works, we also use STAN to generate samples for these parameters using MCMC. The figure below, shows the samples generated by SIR as histograms, and the lines are the samples generated by STAN and MCMC (they agree pretty well). The fourth plot at the bottom right panel shows the data and the regression line.

```## first write the log predictive density function
lpd = function(beta0, beta1, sig){
sigma = exp(sig)
x = lData\$pred
y = lData\$resp
# get the predicted or fitted value
mModMatrix = model.matrix(y ~ x)
mCoef = matrix(c(beta0, beta1), nrow = 2)
iFitted = mModMatrix %*% mCoef # using identity link
## likelihood function with posterior theta
return(sum(dnorm(y, iFitted, sigma, log=T)))
}
```

The scoring function for the model check is the log predictive density, which is basically the log-likelihood function, but using the posterior/fitted parameters.

“In the d-dimensional normal distribution, the logarithm of the density function is a constant plus a $\mathit{X_{d}^{2}}$ distribution divided by -2″ [1],  – and the posterior distribution of the log predictive density has a maximum of  ~ -40.3 and mean of ~ -41.9; with a difference of 1.66, which is close to 3/2 (1.5), predicted from theory, the value of d here is 3 (the number of parameters being estimated). [See Page 171 of Ref 1 for details].

## Predictive Scores

We show three predictive scores along with adjustments for number of parameters being estimated, plus leave one out cross validation. The details for these scores can be found in References [1] and [3], and else where in literature. Lower values of these scores imply a higher predictive accuracy.

### Akaike Information Criterion (AIC)

This is converted to a deviance scale after adjusting for the number of parameters which is 3 in our case.

```> iAIC = (lpd(fit\$mode['beta0'], fit\$mode['beta1'], fit\$mode['sigma']) - 3) * -2
> AIC(fit.lm)
[1] 86.59985
> iAIC
[1] 86.59986
```

The AIC calculated using lpd function and using the AIC function in R on the linear model fit object are the same.

### Deviance Information Criterion (DIC)

This is a somewhat Bayesian version of the AIC, and uses the posterior mean (instead of the maximum likelihood estimate) for θ, and a data-based bias correction.

```> ## pDIC are the effective number of parameters
> ## 2 * [lpd(Expectation(theta)) - Expectation(lpd(Sample of thetas from posterior))]
> # calculate E(lpd(theta))
> eLPD = mean(sapply(1:1000, function(x) lpd(beta0Sample[x], beta1Sample[x], sigSample[x])))
> # calculate lpd(E(theta)) and pDIC
> pDIC = 2 *(lpd(fit\$mode['beta0'], fit\$mode['beta1'], fit\$mode['sigma']) - eLPD)
> iDIC = (lpd(fit\$mode['beta0'], fit\$mode['beta1'], fit\$mode['sigma']) - pDIC) * -2
> pDIC
[1] 3.424087
> iDIC
[1] 87.44804
```

The effective number of parameters will change slightly depending on the simulation size, but both numbers are pretty close (i.e. AIC and DIC).

### Watanabe-Akaike or widely available information criterion (WAIC)

This is a more Bayesian approach and is an approximation to cross-validation. We show one version of this approach here. First calculated log pointwise predictive density, which is slightly different to log predictive density calculated earlier.

```## log pointwise predictive density
lppd = function(beta0, beta1, sig, index){
sigma = exp(sig)
x = lData\$pred[index]
y = lData\$resp[index]
# get the predicted or fitted value
mModMatrix = model.matrix(y ~ x)
mCoef = matrix(c(beta0, beta1), nrow = 2, byrow = T)
iFitted = mModMatrix %*% mCoef # using identity link
## likelihood function with posterior theta
return(mean(dnorm(y, iFitted, sigma, log=F)))
}
# log pointwise predictive probability of the observed data under the fitted
# model is
ilppd = sum(log(sapply(1:15, function(x) lppd(beta0Sample, beta1Sample, sigSample, x))))
```

The above function is slightly different than the lpd function, and works basically on one data point at a time (see in index argument).

```> ## effective numbers of parameters pWAIC1
> pWAIC1 = 2 * (ilppd - eLPD)
>
> iWAIC = -2 * (ilppd - pWAIC1)
> pWAIC1
[1] 2.245018
> iWAIC
[1] 86.26897
```

The effective number of parameters are about 2.2, instead of 3. This is because the effective number of parameters estimated in a model is also a function of the data, and can be considered a random variable – hence in more complex models it is not straight-forward to just count the number of parameters.

### Leave-one-out cross-validation

Cross-validation can be computationally expensive and awkward in structured models, but simply put, the data are partitioned into training and test sets. The parameters are estimated on the training set while the fit is evaluated on the test set.

```## leave one out cross validation
# we need to fit the model 15 times each time removing one data point
lFits = lapply(1:15, function(x){
start = c('sigma' = log(sd(dfData\$VoteShare[-x])), 'beta0'=0, 'beta1'=0)
lData = list(pred=dfData\$IncomeGrowth[-x], resp=dfData\$VoteShare[-x])
laplace(mylogpost, start, lData)
})

# lets take samples for posterior theta
lSamples = lapply(1:15, function(x){
fit = lFits[[x]]
tpar = list(m=fit\$mode, var=fit\$var*2, df=4)
## get a sample directly and using sir (sampling importance resampling with a t proposal density)
lData = list(pred=dfData\$IncomeGrowth[-x], resp=dfData\$VoteShare[-x])
s = sir(mylogpost, tpar, 1000, lData)
return(s)
})

## calculate lppd on the hold out set
iLppd = sapply(1:15, function(x){
s = lSamples[[x]]
log(lppd(s[,'beta0'], s[,'beta1'], s[,'sigma'], x))
})
iLppd = sum(iLppd)
# calculate lppd on all data
# calculated earlier ilppd
# effective number of parameters
iParamCV = ilppd - iLppd
# Given that this model includes two linear coefficients and a variance parameter, these
# all look like reasonable estimates of effective number of parameters. [Gelman 2013]
# on deviance scale iCV is
iCV = -2 * iLppd
# after correction for number of parameters
iCV = -2 * (iLppd - iParamCV)
```

Generally predictive accuracy measures should be used in parallel with posterior predictive checks – starting with simpler models and expanding. This can be applied to comparing nested models: where the full model implements perhaps all the meaningful parameters, and restricted models (where some parameters are restricted or set to 0 or forced to be equal). In this case the complexity of the model and improvement in fit should justify interpretation and additional difficulty in fitting.

## Model Checking: Posterior Predictive Checks

Once a model is fit and parameters estimated, we would look at how well the model explains the data and what aspects of the data generation process in nature are not captured by the model. Most of the material covered in this post follows the examples from:

[1] Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B. (2013). Bayesian Data Analysis, Third Edition (Texts in Statistical Science). Book.

The example data I use here is Simon Newcomb’s experiment to measure the speed of light. The R source code and the data are present in the github repository. Import the data and define a model function i.e. log posterior function to model this data.

```## define a log posterior function
lp = function(theta, data){
# we define the sigma on a log scale as optimizers work better
# if scale parameters are well behaved
s = exp(theta[2])
m = theta[1]
d = data\$vector # observed data vector
log.lik = sum(dnorm(d, m, s, log=T))
log.prior = 1
log.post = log.lik + log.prior
return(log.post)
}

```

Here we assume the data can be modelled using a normal distribution and estimate the parameters.

```## try the laplace function from LearnBayes
fit = laplace(lp, start, lData)

## lets take a sample from this
sigSample.op = rnorm(1000, fit\$mode['sigma'], se['sigma'])
muSample.op = rnorm(1000, mean(lData\$vector), exp(sigSample.op)/sqrt(length(lData\$vector)))

#### the posterior interval reported on P 67 Gelman [2013] is
## y ± 1.997s/ 66 = [23.6, 28.8]
## compare with our intervals
fit\$mode['mu']-1.96*se['mu']; fit\$mode['mu']+1.96*se['mu']
##      mu
## 23.63911
##      mu
## 28.78365
quantile(muSample.op, c(0.025, 0.975))
##     2.5%    97.5%
## 23.61960 28.83061
```

The figure below suggests that the normal model may not be appropriate for this data, as there is an extreme outlier measurement around -44.

There are various approaches to model checking e.g.:

• Using external validation: use data to fit the model (i.e. estimate parameters), make predictions about the new or future data using the model and parameters, collect new real data and compare with the predictions.
• Recycle existing data: if new data is not available. This also involves other modifications to model fitting cross-validation techniques. I will talk about this further in a future post.

How do we compare the predictions with the real data, i.e. we define some metrics to compare against, which are called Test Quantities.

## Prediction of New/Future Data (Posterior Predictive Distribution)

New data or future data can be generated using simulation, from the model and estimated parameters.

```## POSTERIOR PREDICTIVE CHECKING
## observed data should look plausible under the posterior predictive distribution
## Draw simulated values from the joint posterior of Yrep and compare to Yobs and look for systematic differences

## Gelman [2013] P 144 -
## sample 66 values, 20 times, each time drawing a fresh draw of sd and mean from the joint posterior
mDraws = matrix(NA, nrow = 66, ncol=20)

for (i in 1:20){
p = sample(1:1000, size = 1)
s = exp(sigSample.op[p])
m = muSample.op[p]
mDraws[,i] = rnorm(66, m, s)
}

```

In the example shown above we sample the standard deviation parameter, and then sample the mean parameter (conditioned on the standard deviation, notice we use the same index p). We draw 66 samples and repeat the process 20 times, this is akin to repeating the experiment 20 times and each time taking a 66 measurements. The figure below shows some of the histograms of the 20 simulations and the original data.

It appears that the replications do not include the outlier measurement at -44, and suggests that our current model does not capture that part of the data generation process. Using graphical checks can be a bit tedious in high-throughput settings, and we can use some quantitative checks by defining Test Quantities.

## Test Quantities

A test quantity is a function of the original data and replicated data, with some optional additional parameters. We can evaluate the discrepancy in the between the test quantities using the original data and replicated data by calculating a PValue and watch for extreme tail area PValues. We define 5 test quantities in the current case representing: Variance, Mean, Symmetry, Minimum and Maximum.

```# The procedure for carrying out a posterior predictive model check requires specifying a test
# quantity, T (y) or T (y, θ), and an appropriate predictive distribution for the replications
# y rep [Gelman 2008]
## variance
T1_var = function(Y) return(var(Y))

## is the model adequate except for the extreme tails
T1_symmetry = function(Y, th){
Yq = quantile(Y, c(0.90, 0.10))
return(abs(Yq[1]-th) - abs(Yq[2]-th))
}

## min quantity
T1_min = function(Y){
return(min(Y))
}

## max quantity
T1_max = function(Y){
return(max(Y))
}

## mean quantity
T1_mean = function(Y){
return(mean(Y))
}

```
```## calculate bayesian p-value for this test statistic
getPValue = function(Trep, Tobs){
left = sum(Trep <= Tobs)/length(Trep)
right = sum(Trep >= Tobs)/length(Trep)
return(min(left, right))
}

```

Extreme PValues (typically less than 0.01) for the test quantities suggest areas of failure for the model which can be addressed by expanding the model, or ignored if appropriate for the applied question at hand. “The relevant goal is not to answer the question, ‘Do the data come from the assumed model?’ (to which the answer is almost always no), but to quantify the discrepancies between data and model, and assess whether they could have arisen by chance, under the model’s own assumptions.” [1].

Under the normal distribution model, the 5 test quantities show the following PValues.

```> mChecks[,'Normal']
Variance Symmetry      Max      Min     Mean
0.47     0.16     0.01     0.00     0.47
```

The tests for Min and Max quantities fail, suggesting that this model does not perform well in the tail regions. Perhaps using a heavy tailed distribution, like a T distribution with low degrees of freedom or a Contaminated normal distribution will be more useful.

```## define a second log posterior function for mixture with contaminated normal distribution
lp2 = function(theta, data){
# we define the sigma on a log scale as optimizers work better
# if scale parameters are well behaved
s = exp(theta[2])
m = theta[1]
mix = 0.9
cont = theta[3]
d = data\$vector # observed data vector
log.lik = sum(log(dnorm(d, m, s) * mix + dnorm(d, m, s*cont) * (1-mix)))
log.prior = 1
log.post = log.lik + log.prior
return(log.post)
}

# sanity check for function
# choose a starting value
start = c('mu'=mean(ivTime), 'sigma'=log(sd(ivTime)), 'cont'=1)
lp2(start, lData)

## try the laplace function from LearnBayes
fit2 = laplace(lp2, start, lData)
fit2
se2 = sqrt(diag(fit2\$var))

sigSample.op = rnorm(1000, fit2\$mode['sigma'], se2['sigma'])
muSample.op = rnorm(1000, mean(lData\$vector), exp(sigSample.op)/sqrt(length(lData\$vector)))

```

A contaminated normal distribution is a mixture distribution where we add 2 more parameters to the model, the mixing probability and contamination parameter. We fix the mixing probability to 0.9 and track the contamination parameter.

Here are a couple of web links if you are more interested in contaminated normal distributions (Link 1, Link 2).

The third distribution we try is a T distribution with small degrees of freedom. We use a slightly different parameterization of the t-distribution which is nicely explained in another blog here.

```lp3 = function(theta, data){
# function to use to use scale parameter
## see here https://grollchristian.wordpress.com/2013/04/30/students-t-location-scale/
dt_ls = function(x, df, mu, a) 1/a * dt((x - mu)/a, df)
## likelihood function
lf = function(dat, pred){
return(log(dt_ls(dat, nu, pred, sigma)))
}
nu = exp(theta['nu']) ## normality parameter for t distribution
sigma = exp(theta['sigma']) # scale parameter for t distribution
m = theta[1]
d = data\$vector # observed data vector
if (exp(nu) < 1) return(-Inf)
log.lik = sum(lf(d, m))
log.prior = 1
log.post = log.lik + log.prior
return(log.post)
}

# sanity check for function
# choose a starting value
start = c('mu'=mean(ivTime), 'sigma'=log(sd(ivTime)), 'nu'=log(2))
lp3(start, lData)

op = optim(start, lp3, control = list(fnscale = -1), data=lData)
op\$par
exp(op\$par[2:3])

## try the laplace function from LearnBayes
fit3 = laplace(lp3, start, lData)
fit3
se3 = sqrt(diag(fit3\$var))

sigSample.op = rnorm(1000, fit3\$mode['sigma'], se3['sigma'])
nuSample = exp(rnorm(1000, fit3\$mode['nu'], se3['nu']))
# threshold the sample values to above or equal to 1
nuSample[nuSample < 1] = 1

## generate random samples from alternative t-distribution parameterization
## see https://grollchristian.wordpress.com/2013/04/30/students-t-location-scale/
rt_ls <- function(n, df, mu, a) rt(n,df)*a + mu
muSample.op = rnorm(1000, mean(lData\$vector), exp(sigSample.op)/sqrt(length(lData\$vector)))

```

Using the 3 approaches we calculated the test quantities and the PValues for these quantities.

```> round(mChecks, 2)
Normal NormalCont    T
Variance   0.47       0.34 0.34
Symmetry   0.16       0.10 0.14
Max        0.01       0.06 0.18
Min        0.00       0.21 0.14
Mean       0.47       0.48 0.47
```

We can see that the normal model does not seem to fit the data well in terms of outliers, while the Contaminated Normal and T distributions, both fit the data better.

The figure below shows the histogram of the data and the smattering of the density lines from the 200 posterior predictive replications of the data. We can see that the normal model does not fit well in the tails, while the other two show  a reasonable fit, at least it can be explained by sampling variation according to PValues.

While researching the material for this blog I also came across the post about this data set and extreme outliers in data sets at Andrew Gelman’s blog.

## Mixture Distributions and Model Comparison

The following text and code snippets show examples from two books on Bayesian Data Analysis:

[1] Kruschke, J. K. (2014). Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan, second edition. Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan, Second Edition. http://doi.org/10.1016/B978-0-12-405888-0.09999-2

[2] Albert, J., Gentleman, R., Parmigiani, G., & Hornik, K. (2009). Bayesian computation with R. Bayesian Computation with R. http://doi.org/10.1007/978-0-387-92298-0

The examples I refer to here can be seen on Page 50 of Ref 2 (Mixtures of conjugate priors) and Chapter 10 of Ref 1 (Model Comparison .. ). I try and connect the two ideas presented there, where two independent models are either modelled together as a mixture distribution or separately and then compared.

Problem:

We can imagine this problem as two machines or two coins, that have a parameter theta θ, in case of coins the proportion of heads. Coin A is biased to produce more tails, Coin B is biased to produce more heads. We are presented with some data i.e. number of heads and number of tails.

1. Which coin generated this data?
2. What is the distribution of the parameter theta before and after observing the data?

The R code for this analysis and the figures can be found here.

Following the example from ref [2]. Lets define the two distributions:

```### define two models that explain the data
## model m1
## it is a beta distribution with a weight on the left side
g1 = function(theta) dbeta(theta, 6, 14)

## model m2
## it is a beta distribution with a weight on the right side
g2 = function(theta) dbeta(theta, 14, 6)
```

How much weight we assign to each model i.e. the prior probability for each model, and this can be used to generate the mixture distribution.

```## define an array that represents number of models in our parameter space
## each index has a prior weight/probability of being selected
## this can be thought of coming from a categorical distribution
mix.prior = c(50, 50)
mix.prior = mix.prior/sum(mix.prior)

## define a joint prior space for these 2 models
## (mix.prior[1] AND g1) OR (mix.prior[2] AND g2)
g.mix = function(theta) mix.prior[1]*g1(theta) + mix.prior[2]*g2(theta)
```

This function represents the prior distribution or prior beliefs about how the parameter theta for each coin is distributed.

The observed data shows 7 heads and 3 tails. We now define the sampling distribution/likelihood functions and posterior distributions for each model. The posterior distribution represents how much parameter theta changes for each model after observing the data. I have used the word conjugate here, which I will explain at a later time – for the time being think of it as a simple analytical way to calculate the posterior.

```## we flip a coin, perform a binary experiment
## our prior beliefs are that the experiment tends to follow two extreme models
## we will either see more 1s or more 0s, depending on which model (m1 or m2) is
## more representative of the data
data = c(success=7, fail=3)

## the posterior is proportional to Likelihood * Prior
## P(Data | Theta, Model) X P(Theta, Model)

## define the Likelihood function, both models share the same likelihood functional form
lik = function(data, theta) dbinom(data['success'], sum(data), theta)

## since this is a conjugate analysis, as prior is beta distributed and likelihood is binomial
## the posterior can be derived analytically
## using a similar analogy as before when describing priors for each model
## we can define a posterior for each of the 2 models
g1.post = function(theta, data) dbeta(theta, 6+data['success'], 14+data['fail'])
g2.post = function(theta, data) dbeta(theta, 14+data['success'], 6+data['fail'])

```

We have assigned a prior mixture probability in the variable mix.prior, however we need to calculate the posterior mixing probability, which requires some algebra, and I have tried to show that in the code below. For details you have to read the first few pages of Chapter 10 in ref 1.

In order to calculate a mixture probability, i.e. what is the probability each model m1 or m2 after we see the data:

P(Data, Model[1]) = P(Data, Model[1])
P(Model[1] | Data) X P(Data) = P(Data | Model[1]) X P(Model[1])
P(Model[1] | Data) = P(Data | Model[1]) X P(Model[1]) / P(Data) … Equation (1)

where P(Data) can be expanded using summation across all models
P(Data) = Sum for each Model P(Data, Model[i])

We need to calculate a few things:
P(Data | Model[i]) – this is the prior predictive distribution for the data given the selected model. So for each model we solve this equation:

P(Data, Theta) = P(Data | Theta) X P(Theta)

P(Data) = P(Data | Theta) X P(Theta) / P(Theta | Data)
In words this means: Prior predictive distribution for Data = Likelihood X Prior / Posterior

```## Prior predictive probability for Data = Likelihood X Prior / Posterior
## for model 1
data.prior.g1 = function(data, theta){
ret = lik(data, theta) * g1(theta) / g1.post(theta, data)
return(ret)
}
## for model 2
data.prior.g2 = function(data, theta){
ret = lik(data, theta) * g2(theta) / g2.post(theta, data)
return(ret)
}

## P(Data | Model) for each model should be the same for any value of theta
## you can use that as a sanity check
th = seq(0.01, 0.99, by=0.01)
data.g1 = mean(data.prior.g1(data, th))
data.g2 = mean(data.prior.g2(data, th))
```

We have the necessary requirements to solve the equation 1.

```## P(Data) = (P(Data | Model[1]) X P(Model[1])) + (P(Data | Model[2]) X P(Model[2])
## we can calculate the posterior for Model 1
## P(Model[1] | Data) = P(Data | Model[1]) X P(Model[1]) / P(Data)
mix.post = data.g1 * mix.prior[1] / (data.g1 * mix.prior[1] + data.g2 * mix.prior[2])
mix.post = c(mix.post, 1-mix.post)

## (mix.post[1] AND g1.post) OR (mix.post[2] AND g2.post)
g.mix.post = function(theta, data) mix.post[1]*g1.post(theta, data) + mix.post[2]*g2.post(theta, data)

```

The figures below show the posterior distribution of theta for each of the models of the data i.e. head biased and tail biased coins.

You can see that the posterior distribution of theta for model one has shifted more to the right (towards higher values of theta), as the data has influenced this shift. The mixture distribution has also shifted most of its weight towards the right with a very slight bump around 0.4.

We can now approximate this distribution of theta on a grid, and take a random sample from this distribution.

```## approximate the posterior theta on the grid
p.th = g.mix.post(th, data)
p.th = p.th/sum(p.th)
th.sam = sample(th, 10000, replace = T, prob=p.th)
th.sam = th.sam + runif(10000, 0.01/2 * -1, 0.01/2)
summary(th.sam);
##   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
## 0.2041  0.6224  0.6922  0.6748  0.7512  0.9342
```

I will add another post on grid approximations at a later point. For the time being, this is a useful method to generate random samples from a distribution and this sample can be used to calculate various statistics or parameters of the distribution.

We will also approximate this using another technique, where we first define a function that takes one or more parameters as input (defining the parameter space) and the data; and returns one value which is the log posterior at that position in the parameter space. The code below defines two functions, one for each model and uses an optimiser to explore the surface of this function to find the maximum.

```library(LearnBayes)
library(car)
logit.inv = function(p) {exp(p)/(exp(p)+1) }

mylogpost_m1 = function(theta, data){
## theta contains parameters we wish to track
th = logit.inv(theta['theta'])
suc = data['success']
fail = data['fail']

# define likelihood function
lf = function(s, f, t) return(dbinom(s, s+f, t, log=T))

# calculate log posterior
val = lf(suc, fail, th) + dbeta(th, 6, 14, log = T)
return(val)
}

mylogpost_m2 = function(theta, data){
## theta contains parameters we wish to track
th = logit.inv(theta['theta'])
suc = data['success']
fail = data['fail']

# define likelihood function
lf = function(s, f, t) return(dbinom(s, s+f, t, log=T))

# calculate log posterior
val = lf(suc, fail, th) + dbeta(th, 14, 6, log = T)
return(val)
}

# choose sensible starting values for the optimiser
start = c(theta=logit(0.5))
data = c(success=7, fail=3)

fit_m1 = laplace(mylogpost_m1, start, data)
fit_m2 = laplace(mylogpost_m2, start, data)

logit.inv(fit_m1\$mode)
##   theta
##0.428616
logit.inv(fit_m2\$mode)
##   theta
##0.7142694

```

The values for theta representing the maximum of each function are approximately similar to those calculated using grid approximation.

The figure above shows the histogram of the random sample generated using grid approximation, the line represents the value of the mixture posterior g.mix.post, and the red  and green lines represent the maximum points for the 2 models calculated using an optimiser.

## Bayes Factor:

Simply put, the Bayes Factor (BF) indicates how much the prior odds on the two models change after seeing the data. I.e. how much evidence there is for model 1 vs model 2. As usual the technical details can be found in the references 1 and 2.

```## Bayes factor for the ratios of posterior predictive distribution
## of the 2 models
## P(Data | Model[1]) / P(Data | Model[2])
BF = data.g1 / data.g2
## OR posterior odds for the models / prior odds for the 2 models
mix.post[1]/mix.post[2]/(mix.prior[1]/mix.prior[2])
```

The BF calculated analytically is 0.102 in support of model 1. The general convention for a discrete decision about about the models that there is significant evidence for model 1 when BF > 3, and for model 2 when BF < 1/3 (i.e. 0.33). In our case BF for model 2 is 0.1 which shows significant support for model 2, suggesting that the data was generated from model 2 rather than model 1.

The BF can also be approximated using the optimisation approach shown earlier, where the function laplace returns  the P(Data | Model) in the slot \$int.

```> round(exp(fit_m1\$int - fit_m2\$int), 2)
## BF
##   0.09
```

This is very close to the BF of 0.1 calculated analytically. Both analyses suggest that the data was generated by Coin or Machine B.

format: plain


text:
# Form the pair of linear equations in the following problems, and find their solutions graphically.(i) 10 students of class X took part in a Mathematics quiz. If the number of girls is 4 more than the number of boys, find the number of boys and girls who took part in the quiz.(ii) 5 pencils and 7 pens together cost Rs. 50, whereas 7 pencils and 5 pens together cost Rs. 46. Find the cost of one pencil and that of one pen.

#### Complete Python Prime Pack

9 Courses     2 eBooks

#### Artificial Intelligence & Machine Learning Prime Pack

6 Courses     1 eBooks

#### Java Prime Pack

9 Courses     2 eBooks

Given:

10 students of class X took part in the Mathematics quiz. The number of girls is 4 more than the number of boys.

To do:

We have to form a pair of linear equations and solve them graphically. Also, we have to find the number of boys and girls who took part in the quiz.

Solution:

(i) Let the number of girls and boys in the class be $x$ and $y$ respectively.

According to the question,

$x + y = 10$.....(i)

$x - y = 4$.....(ii)

To represent the above equations graphically we need at least two solutions for each of the equations.

For equation $x+y=10$,

$y=10-x$

If $x=5$ then $y=10-5=5$

If $x=7$ then $y=10-7=3$

 $x$ $5$ $7$ $y$ $5$ $3$

For equation $x-y=4$,

$y=x-4$

If $x=4$ then $y=4-4=0$

If $x=6$ then $y=6-4=2$

 $x$ $4$ $6$ $y$ $0$ $2$

The above situation can be plotted graphically as below:

The lines AB and CD represent the equations $x+y=10$ and $x-y=4$ respectively.

The solution of the above equations is the point of intersection.

Therefore,

Number of girls $x=7$

Number of boys $y=3$

The number of boys and girls who took part in the quiz is $3$ and $7$ respectively.

(ii) Let the cost of one pencil and one pen be $x$ and $y$ respectively.

According to the question,

$5x + 7y = 50$.....(i)

$7x + 5y = 46$.....(ii)

To represent the above equations graphically we need at least two solutions for each of the equations.

For equation $5x+7y=50$,

$7y=50-5x$

$y=\frac{50-5x}{7}$

If $x=3$ then $y=\frac{50-5(3)}{7}=\frac{50-15}{7}=\frac{35}{7}=5$

If $x=10$ then $y=\frac{50-5(10)}{7}=\frac{50-50}{7}=\frac{0}{7}=0$

 $x$ $3$ $10$ $y$ $5$ $0$

For equation $7x+5y=46$,

$5y=46-7x$

$y=\frac{46-7x}{5}$

If $x=3$ then $y=\frac{46-7(3)}{5}=\frac{46-21}{5}=\frac{25}{5}=5$

If $x=8$ then $y=\frac{46-7(8)}{5}=\frac{46-56}{5}=\frac{-10}{5}=-2$

 $x$ $3$ $8$ $y$ $5$ $-2$

The above situation can be plotted graphically as below:

The lines AB and CD represent the equations $5x+7y=50$ and $7x+5y=46$ respectively.

The solution of the above equations is the point of intersection.

Therefore,

Cost of one pencil $x=3$

Cost of one pen $y=5$

The cost of one pencil and a pen is $3$ and $5$ respectively.

Updated on 10-Oct-2022 13:19:40


text:
What is the weight per volume method to calculate concentration?

Jan 1, 2015

A weight/volume percent concentration ($\text{w/v%}$) is defined as mass of solute divided by volume of solution and multiplied by 100%. The mathematical expression for a $\text{w/v%}$ percent concentration is

"w/v%" = ("mass of solute")/("volume of solution") * 100%

As an example, a 5% $\text{w/v}$ $N a C l$ solution would have 5 g of $N a C l$ for every 100 mL of solution. A 25% $\text{w/v}$ $N a C l$ solution would have 25 g of $N a C l$ for every 100 mL of solution, and so on.

Weight/volume percent concentrations are characteristic when solids are dissolved in liquids and are often used because volumes are easier to measure than weights.

Another important reason for using such percent concentrations is the fact that dilute solutions have a density that's generally close to 1 g/mL, which makes the volume of the solution expressed in mL nearly numerically equal to the mass of the solution expressed in grams.

As a conclusion, in order to calculate a $\text{w/v}$ percent concentration, you must determine how much solute you have - its mass in grams - and how much volume of solution you have - in mL.


text:
# Thread: Proving that this is a group

1. ## Proving that this is a group

Hello, I would really appreciate any help with the following problem:

Prove that $\displaystyle (G, *)$ is a group, where $\displaystyle G=<-1,1>$ and $\displaystyle a*b=\dfrac{a+b}{1+ab}$,.... $\displaystyle \forall a, b \in G$.

_______________________________

To prove that this is a group, $\displaystyle (G,*)$ must satisfy closure, associativity, existence of an identity element and existence of inverse elements.

Could the property of associativity be demonstrated like this?

$\displaystyle (a*b)*c=(\dfrac{a+b}{1+ab})*c=...=\dfrac{a+b+c+abc }{1+ab+ac+bc}$

$\displaystyle a*(b*c)=a*(\dfrac{b+c}{1+bc})=...=\dfrac{a+b+c+abc }{1+ab+ac+bc}$

Therefore, $\displaystyle (a*b)*c=a*(b*c)$, so the property of associativity is satisfied.

_______________________________

Is the identity element $\displaystyle 0$?

From $\displaystyle a*0= \dfrac{a+0}{1+a \cdot 0}=0*a=a$

it follows $\displaystyle a*0=0*a=a$

_______________________________

Is the inverse element of $\displaystyle a$, $\displaystyle -a$?

Because $\displaystyle a*(-a)=\dfrac{a-a}{1-a^2}=\dfrac{0}{1-a^2}=(a \in <-1,1>)=(-a)*a = 0$

_______________________________

If the above is correct (and I'm not certain that it is), the only property left to be demonstrated is closure, i.e. if $\displaystyle a, b \in <-1,1>$ then $\displaystyle a*b =\dfrac{a+b}{1+ab} \in <-1,1>$.

And here I need your help! How to demonstrate the property of closure?

There is a hint in the text: "Observe that $\displaystyle |1+ab|=1+ab$, and prove that $\displaystyle |a+b| \leq 1+ab$ if and only if $\displaystyle 0 \leq (1-a)(1-b)$ and $\displaystyle 0 \leq (1+a)(1+b)$".

So, how to prove the statements from the hint, and, more importantly, how to apply them to demonstrate the closure property of $\displaystyle (G,*)$?

Many thanks!

2. Originally Posted by gusztav
Hello, I would really appreciate any help with the following problem:

Prove that $\displaystyle (G, *)$ is a group, where $\displaystyle G=<-1,1>$ and $\displaystyle a*b=\dfrac{a+b}{1+ab}$,.... $\displaystyle \forall a, b \in G$.

_______________________________

To prove that this is a group, $\displaystyle (G,*)$ must satisfy closure, associativity, existence of an identity element and existence of inverse elements.

Could the property of associativity be demonstrated like this?

$\displaystyle (a*b)*c=(\dfrac{a+b}{1+ab})*c=...=\dfrac{a+b+c+abc }{1+ab+ac+bc}$

$\displaystyle a*(b*c)=a*(\dfrac{b+c}{1+bc})=...=\dfrac{a+b+c+abc }{1+ab+ac+bc}$

Therefore, $\displaystyle (a*b)*c=a*(b*c)$, so the property of associativity is satisfied.

_______________________________

Is the identity element $\displaystyle 0$?

From $\displaystyle a*0= \dfrac{a+0}{1+a \cdot 0}=0*a=a$

it follows $\displaystyle a*0=0*a=a$

_______________________________

Is the inverse element of $\displaystyle a$, $\displaystyle -a$?

Because $\displaystyle a*(-a)=\dfrac{a-a}{1-a^2}=\dfrac{0}{1-a^2}=(a \in <-1,1>)=(-a)*a = 0$

_______________________________

If the above is correct (and I'm not certain that it is), the only property left to be demonstrated is closure, i.e. if $\displaystyle a, b \in <-1,1>$ then $\displaystyle a*b =\dfrac{a+b}{1+ab} \in <-1,1>$.

And here I need your help! How to demonstrate the property of closure?

There is a hint in the text: "Observe that $\displaystyle |1+ab|=1+ab$, and prove that $\displaystyle |a+b| \leq 1+ab$ if and only if $\displaystyle 0 \leq (1-a)(1-b)$ and $\displaystyle 0 \leq (1+a)(1+b)$".

So, how to prove the statements from the hint, and, more importantly, how to apply them to demonstrate the closure property of $\displaystyle (G,*)$?

Many thanks!

First ,what do you mean by $\displaystyle G=<-1,1>$ ?? The subset of the reals (or complex) with those two elements? If so zero cannot be anything there since it doesn't belong to the set!
Also, the product of $\displaystyle 1\,,\,-1$ isn't defined: $\displaystyle 1*-1:= \frac{1+(-1)}{1+1(-1)}=\frac{0}{0}$ ...!!

If this is so the above isn't a group.

But perhaps you meant the open interval $\displaystyle (-1,1)$ which, by some misterious reason, you denote by $\displaystyle <-1,1>$...(!)

Then you must show $\displaystyle a,b\in(-1,1)\Longrightarrow \frac{a+b}{1+ab}\in (-1,1)\Longleftrightarrow \left|\frac{a+b}{1+ab}\right|<1$ $\displaystyle \Longleftrightarrow (a+b)^2<(1+ab)^2\Longleftrightarrow a^2+2ab+b^2<a^2b^2+2ab+1$

$\displaystyle \Longleftrightarrow a^2(b^2-1)-(b^2-1)>0\Longleftrightarrow (a^2-1)(b^2-1)>0$ , and this last inequality is obviously true.

Tonio

3. Thank you very much, Tonio!

Originally Posted by tonio
But perhaps you meant the open interval $\displaystyle (-1,1)$
Yes, I meant the open interval $\displaystyle \{x \in \mathbb{R} : -1<x<1\}$

Originally Posted by tonio
Then you must show $\displaystyle a,b\in(-1,1)\Longrightarrow \frac{a+b}{1+ab}\in (-1,1)\Longleftrightarrow \left|\frac{a+b}{1+ab}\right|<1$ $\displaystyle \Longleftrightarrow (a+b)^2<(1+ab)^2\Longleftrightarrow a^2+2ab+b^2<a^2b^2+2ab+1$

$\displaystyle \Longleftrightarrow a^2(b^2-1)-(b^2-1)>0\Longleftrightarrow (a^2-1)(b^2-1)>0$ , and this last inequality is obviously true.
Aha! Now everything is clear.
(Because $\displaystyle a \in (-1,1) \Longrightarrow |a|<1 \Longrightarrow a^2<1 \Longrightarrow a^2-1<0$ and, similarly, $\displaystyle b^2-1<0$, so $\displaystyle (a^2-1)( b^2-1)>0$ )


text:
# Polynomials and Factoring

## Adding, Subtracting, and Multiplying Polynomials

To add and subtract polynomials, we "collect like terms", i.e. combine (add or subtract) the coefficients of the terms which have the same variables and exponents.

### Example.

Simplify $3xy+4x^2-5xy+8x^2$.
 $3xy+4x^2-5xy+8x^2$ $=3xy-5xy+4x^2+8x^2$ (Details)Reorganize the order of the terms so those with the same variable expression are next to each other. $=(3-5)xy+(4+8)x^2$ (Details)Combine the coefficients of “like terms” $=-2xy+12x^2$

To multiply two polynomials, we need to take the product of one term in the first polynomial with a term in the second polynomial, repeat this for all such pairs, then add them up.

### Example.

Expand $(x+2)(x^2+5x-1)$. \begin{align*} (x+2)(x^2+5x-1)=& x(x^2)+x(5x)+x(-1)+2(x^2)+2(5x)+2(-1) \\ =&x^3+5x^2-x+2x^2+10x-2 \\ =&x^3+(5+2)x^2+(10-1)x-2 \\ =&x^3+7x^2+9x-2 \end{align*}

Note: Multiplying two terms together involves combining their coefficients and combining their variables. $$\textrm{E.g. }\;\; (2xy^3)(7x^5y^2)=2(7)x(x^5)y^3y^2=14x^6y^5$$ We can also divide polynomials, this is the topic of the next section.


text:
# Finding solution of given recursive equation?

$T(n) = 1+ \sum_{j=0}^{n-1} T(j)$
I've proved $2^n$ be the solution of this equation using induction. But is there any other way to find the solution? I just proved not solved.

• What do you mean by complexity? There seems to be no complexity here, only a function defined by recursion. – Yuval Filmus May 23 '17 at 6:26

The trick in this case is to consider differences: $$T(n)-T(n-1) = \left(1 + \sum_{j=0}^{n-1} T(j)\right) - \left(1 + \sum_{j=0}^{n-2} T(j)\right) = T(n-1).$$ Hence $T(n) = 2T(n-1)$, and so $T(n) = 2^n T(0)$.


text:
# Natural Numbers - Axioms, Addition, Recursion

We see, following the text how to axiomatize the natural numbers and define addition and multiplication from these. We elaborate on only one result, namely one allowing us to make recursive definitions.

Theorem: Let $X$ be a set, let $x_0\in X$ and let $g(n, x)$ be a function $g: \mathbb{N} \times X \to X$. There is a unique function $f: \mathbb{N}\to X$ such that $f(0)= x_0$ and $\forall n \in\mathbb{N}\ f(n++) = g(n, f(n))$.

Proof: The idea is to construct the graph of $f$ as the smallest set containing a specified starting point and closed under a shift.

Let $\mathcal{U} = \{S \subset \mathbb{N}\times X : (0, x_0)\in S,\ (n, x)\in S\Rightarrow (n++, g(n++, x))\in S\}$. Note that $\mathcal{U}$ is non-empty as $\mathbb{N}\times X \in \mathcal{U}$.

Consider the set

It is easy to see that $S_0\in \mathcal{U}$.

Lemma: For every $n\in \mathbb{N}$, there exist a unique $x_n\in X$ such that $(n, x_n) \in S_0$.

Proof: We prove this by induction on $n$. First, let $n=0$. Note that $(0, x_0)\in S$ for all $S\in \mathcal{U}$, and so $(0, x_0) \in S_0$. Now suppose $(0, x_0’)\in S_0$ with $x_0’\neq x_0$. Let $S = \{(n, x) \in S_0: (n, x) \neq (0, x_0’)$. We claim that $S\in \mathcal{U}$.

Firstly, as $(0, x_0) \in S_0$ and $x_0’ \neq x_0$, $(0, x_0)\in S$. Next, suppose $(n, x)\in S \subset S_0$, then $(n++, g(n++, x))\in S_0$. Further, as $n++ \neq 0$, we see that $(n++, g(n++, x)) \neq (0, x_0’)$. Hence $(n++, g(n++, x))\in S$. Thus $S\in \mathcal{U}$.

By definition of $S_0$, as $S\in \mathcal{U}$, $S_0 \subset S$. But $(0, x_0’)$ is in $S_0$ but not $S$, a contradiction. This shows that $x_0\in X$ is the unique point $x \in X$ such that $(0, x)\in S_0$.

The induction step is similar. Assume $(n, x_n)\in S_0$ and this is unique, i.e., $(n, x)\in S_0 \implies x = x_n$. Then $(n, x_n)\in S$ for all $S \in \mathcal{U}$, from which we deduce that $(n++, g(n++, x))\in S_0$. To show uniqueness, suppose $x_{n++}’\in X$ is a point with $(n++, x_{n++}’)\in S_0$ but $x_{n++}’ \neq x_{n++} := g(n++, x)$. Let $S = S_0 \setminus \{(n++, x_{n++}’)\}$. We again show that $S \in \mathcal{U}$ and obtain a contradiction to $S_0 \subset S$.

Namely, as $(0, x_0) \in S_0$ and $n++ \neq 0$, $(0, x_0)\in S$. Next, suppose $(m, x)\in S \subset S_0$, then $(m++, g(m++, x))\in S_0$. If $(m++, g(m++, x)) \notin S$, then we must have $(m++, g(m++, x)) = (n++, x_{n++}’)$. Hence $n = m$, so $x = x_n$. It follows that $x_{n++}’ = g(m++, x) = g(n++, x_n)$, a contradiction. Thus, if $(m, x)\in S$, $(m++, g(m++, x)) \in S$, completing the proof that $S\subset \mathcal{U}$, and contradicting $S_0 \subset S$.

This completes the proof of the lemma.

To prove the theorem, we define the function $f$ by $f(n) = x_n$. As we will see, in terms of set theory the function just corresponds to the set $S_0$. It is straightforward to show by induction than $f$ indeed has the desired properties, and is the unique such function.


text:
# What values of x are not in the domain of the expression (x+5)/(3x^2+23x-36)?

Mar 13, 2017

${D}_{f} \in R - \left\{9 , - \frac{4}{3}\right\}$

#### Explanation:

$\frac{x + 5}{3 {x}^{2} + 23 x - 36}$

The denominator should not be zero.

$3 {x}^{2} + 23 x - 36 \ne 0$

$3 {x}^{2} + 27 x - 4 x - 36 \ne 0$

$\left(x - 9\right) \left(3 x + 4\right) \ne 0$

$x \ne 9$ and $x \ne - \frac{4}{3}$

${D}_{f} \in R - \left\{9 , - \frac{4}{3}\right\}$


text:
# Weierstrass factorization theorem

Weierstrass factorization theorem

In mathematics, the Weierstrass factorization theorem in complex analysis, named after Karl Weierstrass, asserts that entire functions can be represented by a product involving their zeroes. In addition, every sequence tending to infinity has an associated entire function with zeroes at precisely the points of that sequence.

A second form extended to meromorphic functions allows one to consider a given meromorphic function as a product of three factors: the function's poles, zeroes, and an associated non-zero holomorphic function.

Motivation

The consequences of the fundamental theorem of algebra are twofold.citation|last=Knopp|first=K.|contribution=Weierstrass's Factor-Theorem|title=Theory of Functions, Part II|location=New York|publisher=Dover|pages=1–7|year=1996.] Firstly, any finite sequence,$\left\{c_n\right\}$, in the complex plane has an associated polynomial that has zeroes precisely at the points of that sequence:

:$,prod_n \left(z-c_n\right).$

Secondly, any polynomial function in the complex plane, $p\left(z\right)$, has a factorization

:$,p\left(z\right)=aprod_n\left(z-c_n\right),$

where "a" is a non-zero constant and "c""n" are the zeroes of "p".

The two forms of the Weierstrass factorization theorem can be thought of as extensions of the above to entire functions. The necessity of extra machinery is demonstrated when one considers whether the product

:$,prod_n \left(z-c_n\right)$

defines an entire function if the sequence, $\left\{c_n\right\}$, is not finite. The answer is never, because the now-infinite product will not converge. Thus one cannot, in general, define an entire function from a sequence of prescribed zeroes or represent an entire function by its zeroes using the expressions yielded by the fundamental theorem of algebra.

A necessary condition for convergence of the infinite product in question is: each factor, $\left(z-c_n\right)$, must approach 1 as $n oinfty$. So it stands to reason that one should seek a function that could be 0 at a prescribed point, yet remain near 1 when not at that point and furthermore introduce no more zeroes than those prescribed. Enter the genius of Weierstrass' "elementary factors". These factors serve the same purpose as the factors, $\left(z-c_n\right)$, above.

The elementary factors

These are also referred to as "primary factors".citation|last=Boas|first=R. P.|title=Entire Functions|publisher=Academic Press Inc.|location=New York|year=1954|isbn=0821845055|oclc=6487790, chapter 2.]

For $n in mathbb\left\{N\right\}$, define the "elementary factors":citation|last=Rudin|first=W.|title=Real and Complex Analysis|edition=3rd|publisher=McGraw Hill|location=Boston|pages=301–304|year=1987|isbn=0070542341|oclc=13093736.]

:

Their utility lies in the following lemma:

Lemma (15.8, Rudin) for |"z"| ≤ 1, "n" ∈ No

:$vert 1 - E_n\left(z\right) vert leq vert z vert^\left\{n+1\right\}.$

The two forms of the theorem

equences define holomorphic functions

Sometimes called the Weierstrass theoremMathWorld | urlname=WeierstrasssTheorem | title=Weierstrass's Theorem]

If $lbrace z_i brace_i subset mathbb\left\{C\right\}-\left\{0\right\}$ is a sequence such that:
#$vert z_i vert ightarrow infty$ as $i ightarrow infty$
#there is a sequence, $lbrace p_i brace_i subset mathbb\left\{N\right\}_o$, such that for all "r" > 0, $sum_\left\{i\right\} left\left( frac\left\{r\right\}\left\{vert z_i vert\right\} ight\right)^\left\{1+p_i\right\} < infty.$Then there exists an entire function that has (only) zeroes at every point of $lbrace z_i brace$; in particular, "P" is such a function:

: $P\left(z\right)=prod_\left\{i=1\right\}^infty E_\left\{p_i\right\}left\left(frac\left\{z\right\}\left\{z_i\right\} ight\right).$

*The theorem generalizes to: sequences in open subsets (and hence regions) of the Riemann sphere have associated functions that are holomorphic in those subsets and have zeroes at the points of the sequence.
*Note also that the case given by the fundamental theorem of algebra is incorporated here. If the sequence, $\left\{ z_i \right\}$ is finite then setting $p_i = 0$ suffices for convergence in condition 2, and we obtain: $, P\left(z\right) = prod_n \left(z-z_n\right)$.

Holomorphic functions can be factored

Sometimes called the Weierstrass Product/Factor/Factorization theorem.MathWorld | urlname=WeierstrassProductTheorem | title=Weierstrass Product Theorem] Sometimes called the Hadamard Factorization theorem; for example c.f. Boas.

If "f" is a function holomorphic in a region, $Omega$, with zeroes at every point of $lbrace z_i brace_i subset mathbb\left\{C\right\}-\left\{0\right\}$then there exists an entire function "g", and a sequence $lbrace p_i brace_i subset mathbb\left\{R\right\}_o^+$ such that:

$f\left(z\right)=e^\left\{g\left(z\right)\right\} prod_\left\{i=1\right\}^infty E_\left\{p_i\right\}left\left(frac\left\{z\right\}\left\{z_i\right\} ight\right)$

References

* Entire function

Wikimedia Foundation. 2010.

### Look at other dictionaries:

• Weierstrass preparation theorem — In mathematics, the Weierstrass preparation theorem is a tool for dealing with analytic functions of several complex variables, at a given point P. It states that such a function is, up to multiplication by a function not zero at P, a polynomial… …   Wikipedia

• Weierstrass theorem — Several theorems are named after Karl Weierstrass. These include: *The Weierstrass approximation theorem, also known as the Stone Weierstrauss theorem *The Bolzano Weierstrass theorem, which ensures compactness of closed and bounded sets in R n… …   Wikipedia

• Karl Weierstrass — Infobox Scientist name = Karl Weierstrass |300px caption = Karl Theodor Wilhelm Weierstrass (Weierstraß) birth date = birth date|1815|10|31|mf=y birth place = Ostenfelde, Westphalia death date = death date and age|1897|2|19|1815|10|31|mf=y death… …   Wikipedia

• Mittag-Leffler's theorem — In complex analysis, Mittag Leffler s theorem concerns the existence of meromorphic functions with prescribed poles. It is sister to the Weierstrass factorization theorem, which asserts existence of holomorphic functions with prescribed zeros. It …   Wikipedia

• List of mathematics articles (W) — NOTOC Wad Wadge hierarchy Wagstaff prime Wald test Wald Wolfowitz runs test Wald s equation Waldhausen category Wall Sun Sun prime Wallenius noncentral hypergeometric distribution Wallis product Wallman compactification Wallpaper group Walrasian… …   Wikipedia

• Gamma function — For the gamma function of ordinals, see Veblen function. The gamma function along part of the real axis In mathematics, the gamma function (represented by the capital Greek letter Γ) is an extension of the factorial function, with its… …   Wikipedia

• List of complex analysis topics — Complex analysis, traditionally known as the theory of functions of a complex variable, is the branch of mathematics that investigates functions of complex numbers. It is useful in many branches of mathematics, including number theory and applied …   Wikipedia

• Entire function — In complex analysis, an entire function, also called an integral function, is a complex valued function that is holomorphic over the whole complex plane. Typical examples of entire functions are the polynomials and the exponential function, and… …   Wikipedia

• Zero (complex analysis) — In complex analysis, a zero of a holomorphic function f is a complex number a such that f(a) = 0. Contents 1 Multiplicity of a zero 2 Existence of zeros 3 Properties …   Wikipedia

• Polygamma function — In mathematics, the polygamma function of order m is defined as the ( m + 1)th derivative of the logarithm of the gamma function::psi^{(m)}(z) = left(frac{d}{dz} ight)^m psi(z) = left(frac{d}{dz} ight)^{m+1} lnGamma(z).Here :psi(z) =psi^{(0)}(z) …   Wikipedia


text:
# 9 Model definition

IQR Tools implements a general and powerful ODE and biochemical reaction expression based model description language, allowing to describe most models used in systems biology, systems pharmacology, and pharmacometrics. Tools are available to support the development, manipulation, and simulation of models.

## 9.1 Model definition basics

### 9.1.1 Biochemical reaction model

The basics of model representation in IQR Tools is illustrated for a model describing the (irreversible) formation of a complex C between a substrate A and an enzyme B (see Figure below). These types of processes are common as motifs in larger models of signaling pathways.

In IQR Tools models are described in txt-files containing different sections for the different elements (e.g., states, parameters, etc.) of the model. Some of the sections may be empty. MODEL FUNCTIONS and MODEL EVENTS for example are used for representation of more complex models which is explained in Section Advanced model representation. The txt-file of the model takes the following form:

********** MODEL NAME
Biochemical reaction model
********** MODEL NOTES
Concentration changes due to biochemical reaction given as differential equations
********** MODEL STATES
d/dt(Substrate) = -R1
d/dt(Enzyme)    = -R1
d/dt(Complex)   =  R1

Substrate(0) = 1
Enzyme(0)    = 1
Complex(0)   = 0
********** MODEL PARAMETERS
k1 = 0.5
********** MODEL VARIABLES

********** MODEL REACTIONS
R1 = k1 * Substrate * Enzyme
********** MODEL FUNCTIONS

********** MODEL EVENTS

The MODEL STATES section contains the ordinary differential equations that define the changes of each state with time and their initial value at t=0. Initial state values may be omitted if they are 0. Another feature is that an initial state does not need to be a numerical value but could also be a model parameter.

The reaction rate for this binding reaction is defined in the MODEL REACTIONS section. It could be defined directly in the state equations, but this may be confusing for more complex reaction networks.

Having defined the model in a text file, it can be loaded, simulated, and the simulation results visualized with IQR Tools. The chapter Basic simulation of models illustrates how to simulate models with IQR tools in more detail.

# load model
filename <- "material/01-02-ModelImplementation/model_Example1.txt"
model <- IQRmodel(filename)

# simulate the model for 100 time units
sim_results <- sim_IQRmodel(model,100)

# plot the simulation results for the complex
plot(sim_results)

Alternatively to the differential equations, a representation can be used that resembles the common syntax for writing (bio-)chemical reactions.

********** MODEL NAME
Biochemical reaction model

********** MODEL NOTES
Biochemical reaction representation
********** MODEL STATE INFORMATION
Substrate(0) = 1
Enzyme(0)    = 1
Complex(0)   = 0
********** MODEL PARAMETERS
k1 = 0.5
********** MODEL VARIABLES

********** MODEL REACTIONS
Substrate + Enzyme  => Complex  :  R1
vf = kon * Drug * Target
********** MODEL FUNCTIONS
********** MODEL EVENTS

Here, the MODEL STATES section was substituted by the MODEL STATE INFORMATION section, containing only the initial state values. In the MODEL REACTIONS section, the reaction equation is given followed by its name (R1) after a colon. The (forward) rate (vf) is defined below.

In this representation it is easy to switch the reaction to be reversible by adjusting the arrow type and giving the rate for the reverse reaction:

********** MODEL REACTIONS
Substrate + Enzyme  <=> Complex  :  R1
vf = kon * Drug * Target
vr = koff * Complex

### 9.1.2 One compartment linear PK

A one-compartment linear pharmacokinetic model with an initial single dose can be implented in IQR Tools as follows:

********** MODEL NAME
Linear 1-compartmental PK model
********** MODEL NOTES
Linear 1-compartmental PK model with initial single dose
********** MODEL STATES
d/dt(Ac) = -CL/Vc*Ac
Ac(0) = Dose
********** MODEL PARAMETERS
CL    =  20    # (L/hour)   Clearance
Vc    =  32    # (L)         Central volume
Dose  = 100    # (mg)       Dose
********** MODEL VARIABLES
Cc = Ac/Vc * 1000     # (ng/mL) Plasma concentration
********** MODEL REACTIONS
********** MODEL FUNCTIONS
********** MODEL EVENTS

The MODEL VARIABLES section is here used to define the plasma concentration based on the amount (a state) and the central volume (a parameter). Variables could also be used in the definiton of the ODEs in the MODEL STATES section or in the MODEL REACTIONS section.

# load model
filename <- "material/01-02-ModelImplementation/model_Example2.txt"
model <- IQRmodel(filename)

# simulate model for 20 hours
sim_results <- sim_IQRmodel(model,20)

# plot the simulated plasma concentration
plot(sim_results)

### 9.1.3 PK/PD

ODEs and biochemical reactions can be used simultaneously. As an example, the pharmacokinetic model is linked with a binding reaction to a clinical response.

The PK and the changes of the clinical symptoms in dependency of a target molecule are coded with differential equations whereas the binding reaction of the drug to the target molecule is defined as a biochemical reaction.

********** MODEL NAME
Drug target binding with indirect clinical response
********** MODEL NOTES
Target indirectly stimulating clinical symptoms.
Drug binds to target forming a non-stimulating complex.
********** MODEL STATE INFORMATION
d/dt(Ac)            = -CL/Vc*Ac
d/dt(ClinSym)       = rinS - koutS * ClinSym
Ac(0)               =  Dose
ClinSym(0)          =  ClinSymBase
Target(0)           =  TargetBase
********** MODEL PARAMETERS
# PK
CL          =  20       # (L/hour)   Clearance
Vc          =  32       # (L)        Central volume
Dose        =  10       # (mg)       Dose
MW          = 400       # (g/mol)    Molecular weight
TargetBase  =   3       # (nM)       Baseline target levels
ClinSymBase = 100       # (score)    Clinical symptom score at baseline
koutS       =   0.2     # (-)        Turnover rate of clinical symptoms
KD          =  50       # (nM)       Binding dissociation constant
koff        =   0.1     # (1/h)      Dissociation rate constant
********** MODEL VARIABLES
Cc        = Ac/Vc * 1000                     # (ng/mL) plasma concentration
Stim      = ClinSymBase/TargetBase * koutS   # clinical symptom stimulation factor
rinS      = Stim * Target                    # symptom input rate
kon       = koff/KD                          # binding rate constant
********** MODEL REACTIONS
# Reversible binding reaction
Cc + Target  <=> Complex  :  Rbind
vf = kon * Cc/MW*1000 * Target
vr = koff * Complex
********** MODEL FUNCTIONS
********** MODEL EVENTS

The PK, free target molecule levels and clinical symptoms are simulated in IQR Tools:

## 9.2 Dosing representation

Usually, pharmacokinetic models need to handle different and complex dosing regimens. An implementation directly in the model as shown in the PK model example above is not sufficient for more complex regimens or other tasks than simulation. Instead, the PK model and the dosing are separated. In the model, a variable INPUT1 is added to the compartment (amount of drug) to which the dose is administered.

********** MODEL NAME
Linear 1-compartmental PK model
********** MODEL NOTES
Linear 1-compartmental PK model with initial single dose
********** MODEL STATES
d/dt(Ac) = -CL/Vc*Ac    + INPUT1
Ac(0) = 0
********** MODEL PARAMETERS
CL    =  20    # (L/hour)   Clearance
Vc    =  32    # (L)         Central volume
********** MODEL VARIABLES
Cc = Ac/Vc * 1000     # (ng/mL) Plasma concentration
********** MODEL REACTIONS
********** MODEL FUNCTIONS
********** MODEL EVENTS

The dosing is defined by the function IQRdosing() with (at least) defining the time of dosing (“TIME”), the dose amount (“AMT”) and the input number (“ADM”) as there can be more input variables present. This way it is easy to simulate the PK model with two bolus inputs of 10 mg at time 0 and 10.

# load model
filename <- "material/01-02-ModelImplementation/model_Example4.txt"
model <- IQRmodel(filename)

# define the dosing scheme

# simulate model for 20 time units
sim_results <- sim_IQRmodel(model,20,dosingTable=dosing)

# plot the simulation results for the central compartment
plot(sim_results)

Infusions (and zero-order absorption after oral dosing) are implemented by defing a rate (RATE) of infusion or the infusion duration (TINF). First order absorption would be coded in the model with a bolus administration to an additional dosing compartment. For any dosing event, the parameter “LAGTIME” can specify a lag time before the dosing.

### 9.3.1 Example

Assume a drug which basically has a 1-compartmental PK as in the example above, but a limited clearance and parameters changing with time and impacted by the dosing. A PD effect (PD) is directly linked to the plasma concentration (c).

To implement this model, (advanced) model syntax is used. Conditional statements are defined for modulation of the absorption rate constant. Model events are used to define when the last dose was given. Interpolation function is used to define the time-varying Michaelis-Menten constant. The Michealis-Menten equation is a user-defined function.

An implementation of the model as an IQR model is:

********** MODEL NAME
********** MODEL NOTES
1-compartmentl PK model with ka doubled in the first 30 min after dosing and a Michaelis-Menten like limitation on clearance that increases with time.
Assume time in hours.
********** MODEL STATES
# PK model
# dummy state indicating when to apply increased KA
d/dt(DUMMY1) = INPUT1
# dummy state indicating whether first dose was already given
d/dt(DUMMY2) = INPUT1
# count time after dose only when first dose was already given
d/dt(timeAfterDose) = flag
Ac(0) = 0
********** MODEL PARAMETERS
ka = 0.6
CL = 20
Vc = 12
********** MODEL VARIABLES
# Set instantaneous KA based on dummy state 1
#  (DUMMY1 is not empty)
kaINST = piecewise(ka,eq(DUMMY1,0),2*ka)

# Set flag if first dose was already given
#  (DUMMY2 is not empty)
flag = piecewise(0,eq(DUMMY2,0),1)

# Time varying Michaelis-Menten constant
kM = interp1([0,12,24],[0.1,1,10],time)

# Plasma concentration
Cc = Ac/Vc

# PD effects
PD = log(Cc)^2
********** MODEL REACTIONS
CLinst = CL/Vc * MM_function(Cc,kM)
********** MODEL FUNCTIONS
# Define Michaelis-Menten function
MM_function(x,kA) = x/(x+kA)
********** MODEL EVENTS
# Reset time after dose to 0 when dosing occurs
Event1 = gt(DUMMY1,0),timeAfterDose,0

# Empty dummy state after 0.5 hours
Event2 = gt(timeAfterDose,0.5),DUMMY1,0

These are the simulated timecourses for 5 doses of 10 mg each 12 hours:

# load model
filename <- "material/01-02-ModelImplementation/model_Example5.txt"
model <- IQRmodel(filename)

# define the dosing scheme

# simulate model for 20 time units
sim_results <- sim_IQRmodel(model,72,dosingTable=dosing)

# plot the simulation results for the central compartment
plot(sim_results)

In the following, the different parts of the model description above are discussed in detail. For general information about the parts of the model description, please refer to the model section descripton in the technical details.

### 9.3.2 Lagtimes

Whenever an INPUT* is defined in a model (where “*” is the number of the input) then a parameter Tlag*, which represents the lagtime for this input, is assumed to be present in the model. The user can code this parameter manually in the text representation of the model, or if this parameter is not present, then it is added automatically upon import of a text representation, and initialized with 0.

### 9.3.3 Mathematical functions

Basic mathematical functions, i.e., the natural logarithm (“log”) and power (“^”), are needed for linking the the plasma concentrations to the PD effect. Many other basic mathematical functions are available with an intuitive notation. Most commonly used are

• log (natural logarithms)
• exp (natural exponential function)
• min (minimum)
• max (maximum)
• sqrt (square root)
• ^ (power operator)

Please refer to the section Pre-defined functions for a complete list.

### 9.3.4 Implementing conditional statements (if-then-else)

In the model in section Advanced model representation the absorption rate value used in the ODE is high if the last dose was given less than 30 min ago and low otherwise. This is a conditional statements (if-then-else) that is implemented in IQRtool via the “piecewise” function:

# Set instantaneous KA based on dummy state 1
#  (DUMMY1 is not empty)
kaINST = piecewise(ka,eq(DUMMY1,0),2*ka)

Here, DUMMY1 indicates the time after dose until 30 min (DUMMY1 is 1 for 0-30 min after the last dose and 0 otherwise). The second argument checks whether DUMMY1 is 0. If this is true, the first argument (ka) is returned, the third argument (2*ka) otherwise.

Conditional statements use value comparisons (e.g., “a is less than b”, “c is equal to b”) and Boolean operators (AND, OR, etc. ) as also shown in the example above which are implemented using the IQR Tools function in the following table:

Comparison/Operator IQR Tools function
< lt
> gt
>= le
>= ge
= eq
AND and
OR or

### 9.3.5 Interpolation functions

In the model in section Advanced model representation the Michaelis-Menten constant is linearily interpolated from 0.1 to 1 and 1 to 10 between time points 0 and 12 hours, and 12 and 24 hours. The first argument to the interpolation function “interp1” defines the time points and the second which values kM should assume at these timepoints. The third argument indicates the variable specifying the current time point.

# Time varying Michaelis-Menton constant
kM = interp1([0,12,24],[0.1,1,10],time)

The table below lists the interpolation functions available in IQR tools:

Function Description
interp0 Zero order interpolation (constant values between support points)
interp1 First order (linear) interpolation
interpcs Cubic spline interpolation

### 9.3.6 MODEL FUNCTIONS section

The user can also define own functions that can be used in the other model sections. In the example in section Advanced model representation, the Michaelis-Menten equation was defined as a function such that could be called by its name in the MODEL REACTIONS section:

********** MODEL REACTIONS
CLinst = CL/Vc * MM_function(Cc,kM)

********** MODEL FUNCTIONS
# Define Michaelis-Menten function
MM_function(x,kA) = x/(x+kA)

Definition of functions can be useful for improving readability in case of complicated expressions or terms that are repeatedly used in the model.

### 9.3.7 MODEL EVENTS section

In the model in section Advanced model representation, the state “DUMMY1” is an indicator of whether the last dosing time point was no more than 30 min ago or not. It is recognizing a dose as the INPUT1 is added to this state and DUMMY1 is no longer empty after a dose. However, it is reset to 0 at 0.5h so that it can be changed also for later doses. At the same time, the state “timeafterdose” which is increasing along with time needs to be reset to zero if a dose is given, i.e., when DUMMY1 is not empty anymore.

In the IQR model representation, instant changes of states (or variables) are implemented using model events:

********** MODEL EVENTS
# Reset time after dose to 0 when dosing occurs
Event1 = gt(DUMMY1,0),timeAfterDose,0

# Empty dummy state after 0.5 hours
Event2 = gt(timeAfterDose,0.5),DUMMY1,0

An event statement consists of a condition when the event occurs followed by the name of the state to be changed and the value to which it is reset. The “Event1” in the example reads: if “DUMMY1” becomes greater than 0, set “timeafterdose” to 0.

## 9.4 Handling of models in R

### 9.4.1 Model import

An example of how to import an existing model using the function IQRmodel() is shown below. The extension of the model file name is arbitrary. It is suggested to use “.txt” for ODE based models, and “.txtbc” for models in the biochemical representation. But other extensions can be used if desired.

When importing a model an object of class “IQRmodel” is created.

# Load existing model
filename <- "material/01-02-ModelImplementation/model_Example1.txt"
model <- IQRmodel(filename)

# Inspect class of created object
class(model)
## [1] "IQRmodel"

### 9.4.2 Support of SBML

Systems Pharmacology gains in interest in the Modeling and Simulation community. Thus does also the Systems Biology Markup Language (SBML), in which many systems pharmacology type of models are coded (see, e.g. Biomodels Database).

SBML, however, is very badly supported when using R. Some packages (such as SBMLR) are available and claim to be able to import and simulate SBML models. However, typically, not all features of SBML models are considered.

In order to allow IQR Tools to import SBML models a workaround needed to be found. This workaround comes in the form of the IQRsbml software that needs to be installed separately on the computer.

Once installed, the path to the executable needs to be provided in the setup of IQR Tools (see Section ??. After that IQR Tools can easily import SBML models (Level 2 Version 1,2,3,4) as IQRmodels with the following command:

modelSBML <- IQRmodel("sbmlModelFile.xml")

IQRsbml is only provided for Windows 64 bit systems.

### 9.4.3 Basic model information

Basic information about the model is displayed when calling the model variable. The IQRmodel object stores the definitions given in the model sections and some further elements as list elements that can be accessed using the $-operator. Further information is stored as attributes which are described in the model attributes section in the model definition chapter. # Display basic model information model ## IQRmodel ## ======== ## Name: Biochemical reaction model ## Number States: 3 ## Number Parameters: 1 ## Number Variables: 0 ## Number Reactions: 1 ## Number Functions: 0 ## Number Events: 0 ## Number Inputs: 0 ## Number Outputs: 0 # Names of model sections/elements names(model) ## [1] "name" "notes" "parameters" "states" "reactions" # Access a model section, e.g., the model parameters model$parameters
## $k1 ##$k1$value ## [1] 0.5 ## ##$k1$type ## NULL ## ##$k1$compartment ## NULL ## ##$k1$unittype ## NULL ## ##$k1\$notes
## NULL

### 9.4.4 Model export

To export a model to a txt file we use the export_IQRmodel() function. In the following example the model is saved in txt file called “model.txt” in the current working directory. The model can also be exported in the biochemical reaction representation by setting “FLAGbc = TRUE”. The export function can also be used to view the model on the console during working.

# export an existing model
filename <- file.path(getwd(),"model.txt")
export_IQRmodel(model,filename)
# View the exported model on the console
cat(export_IQRmodel(model))
## ********** MODEL NAME
##
## Biochemical reaction model
##
## ********** MODEL NOTES
##
## Concentration changes due to biochemical reaction given as differential equations
##
## ********** MODEL STATES
##
## d/dt(Substrate) = -R1
## d/dt(Enzyme) = -R1
## d/dt(Complex) = R1
##
## Substrate(0) = 1
## Enzyme(0) = 1
## Complex(0) = 0
##
## ********** MODEL PARAMETERS
##
## k1 = 0.5
##
## ********** MODEL VARIABLES
##
##
## ********** MODEL REACTIONS
##
## R1 = k1 * Substrate * Enzyme
##
## ********** MODEL FUNCTIONS
##
##
## ********** MODEL EVENTS
##
##
## ********** MODEL INITIAL ASSIGNMENTS
##
##
# View the exported model as biochemical reactions on the console
cat(export_IQRmodel(model, FLAGbc = TRUE))
## For export in BC notation the underlying assumptions are:
##   * All reaction rates are defined in amount/time units
##   * Species ideally in amount and ODEs defined by sum of reaction rates
##     * Stoichiometric factors can be numeric or parameters and need to be define in front of reaction names
##       Example: d/dt(speciesAmount) = 2*reaction1 - reaction2 - 1*reaction3
## ********** MODEL NAME
##
## Biochemical reaction model
##
## ********** MODEL NOTES
##
## Concentration changes due to biochemical reaction given as differential equations
##
## ********** MODEL STATE INFORMATION
##
##
## Substrate(0) = 1
## Enzyme(0) = 1
## Complex(0) = 0
##
## ********** MODEL PARAMETERS
##
## k1 = 0.5
##
## ********** MODEL VARIABLES
##
##
## ********** MODEL REACTIONS
##
## Substrate+Enzyme => Complex : R1
##  vf = k1 * Substrate * Enzyme
##
## ********** MODEL FUNCTIONS
##
##
## ********** MODEL EVENTS
##
##
## ********** MODEL INITIAL ASSIGNMENTS
##
##

## 9.5 PK model library

A PK model library is available in the documentation material: material/ModelLibraryPK.

## 9.6 Example models

In this section some model definitions for certain models are shown.

### 9.6.1 PBPK

********** MODEL NAME
Example PBPK model
********** MODEL NOTES
Dose: mg
Concentrations: ng/mL
Time: hours
********** MODEL STATES
# Differential equations, implementing a PBPK model
d/dt(Lungs)          = QLU*(Venous_Blood/VVB - Lungs/KbLU/VLU)
d/dt(Heart)          = QHT*(Arterial_Blood/VAB - Heart/KbHT/VHT)
d/dt(Brain)          = QBR*(Arterial_Blood/VAB - Brain/KbBR/VBR)
d/dt(Muscles)        = QMU*(Arterial_Blood/VAB - Muscles/KbMU/VMU)
d/dt(Skin)           = QSK*(Arterial_Blood/VAB - Skin/KbSK/VSK)
d/dt(Spleen)         = QSP*(Arterial_Blood/VAB - Spleen/KbSP/VSP)
d/dt(Pancreas)       = QPA*(Arterial_Blood/VAB - Pancreas/KbPA/VPA)
d/dt(Liver)          = QHA*Arterial_Blood/VAB + QSP*Spleen/KbSP/VSP + QPA*Pancreas/KbPA/VPA + QST*Stomach/KbST/VST + QGU*Gut/KbGU/VGU - CLint*fub*Liver/KbLI/VLI - QLI*Liver/KbLI/VLI
d/dt(Stomach)        = QST*(Arterial_Blood/VAB - Stomach/KbST/VST)
d/dt(Gut)            = QGU*(Arterial_Blood/VAB - Gut/KbGU/VGU)
d/dt(Bones)          = QBO*(Arterial_Blood/VAB - Bones/KbBO/VBO)
d/dt(Kidneys)        = QKI*(Arterial_Blood/VAB - Kidneys/KbKI/VKI)
d/dt(Arterial_Blood) = QLU*(Lungs/KbLU/VLU - Arterial_Blood/VAB)
d/dt(Venous_Blood)   = QHT*Heart/KbHT/VHT + QBR*Brain/KbBR/VBR + QMU*Muscles/KbMU/VMU + QAD*Adipose/KbAD/VAD + QSK*Skin/KbSK/VSK + QLI*Liver/KbLI/VLI + QBO*Bones/KbBO/VBO + QKI*Kidneys/KbKI/VKI + QRB*Rest_of_Body/KbRB/VRB - QLU*Venous_Blood/VVB + INPUT1
d/dt(Rest_of_Body)   = QRB*(Arterial_Blood/VAB - Rest_of_Body/KbRB/VRB)

********** MODEL PARAMETERS
# Define some parameters
# These are the ones that are estimated in the example
CLint   = 1998.196  # Intrinsic clearance
KbBR    = 3.004166  # Partition coefficient - Brain
KbMU    = 1.349859  # Partition coefficient - Muscle
KbBO    = 1.030455  # Partition coefficient - Bone
KbRB    = 1.349859  # Partition coefficient - Rest of body

BP      = 0.61     # Blood:plasma partition coefficient
fup     = 0.028    # Fraction unbound in plasma

KbLU    = 2.301129
KbHT    = 3.066387
KbSK    = 0.5922657
KbSP    = 1.380437
KbPA    = 1.380437
KbLI    = 5.814763
KbST    = 1.380437
KbGU    = 3.32876
KbKI    = 3.732581

# Define regression parameters
WT      = 70
********** MODEL VARIABLES
# Determine fraction unbound in blood
fub     = fup/BP

# Define regional blood flows
CO      = (187.00*WT^0.81)*60/1000 # Cardiac output (L/h) from White et al (1968)
QHT     = 4.0 *CO/100
QBR     = 12.0*CO/100
QMU     = 17.0*CO/100
QSK     = 5.0 *CO/100
QSP     = 3.0 *CO/100
QPA     = 1.0 *CO/100
QLI     = 25.5*CO/100
QST     = 1.0 *CO/100
QGU     = 14.0*CO/100
QHA     = QLI - (QSP + QPA + QST + QGU) # Hepatic artery blood flow
QBO     = 5.0 *CO/100
QKI     = 19.0*CO/100
QRB     = CO - (QHT + QBR + QMU + QAD + QSK + QLI + QBO + QKI)
QLU     = QHT + QBR + QMU + QAD + QSK + QLI + QBO + QKI + QRB

# Define organ volume
# Organ volume = organ weight / organ density
VLU     = (0.76 *WT/100)/1.051
VHT     = (0.47 *WT/100)/1.030
VBR     = (2.00 *WT/100)/1.036
VMU     = (40.00*WT/100)/1.041
VSK     = (3.71 *WT/100)/1.116
VSP     = (0.26 *WT/100)/1.054
VPA     = (0.14 *WT/100)/1.045
VLI     = (2.57 *WT/100)/1.040
VST     = (0.21 *WT/100)/1.050
VGU     = (1.44 *WT/100)/1.043
VBO     = (14.29*WT/100)/1.990
VKI     = (0.44 *WT/100)/1.050
VAB     = (2.81 *WT/100)/1.040
VVB     = (5.62 *WT/100)/1.040
VRB     = (3.86 *WT/100)/1.040

# Determine venous blood concentration in ug/mL
CVBngmL = 1000 * Venous_Blood / (VVB*BP)

# Define output for estimation
tOut    = log(CVBngmL)

# Output for fitting
OUTPUT1 = tOut # Log of venous blood concentration [ug/mL]
********** MODEL REACTIONS
********** MODEL FUNCTIONS
********** MODEL EVENTS

### 9.6.2 Friberg neutropenia

********** MODEL NAME
Friberg neutropenia
********** MODEL NOTES
Example neutropenia model
Dose: mg
Output unit: 10^9 cells/L
Time: hours
********** MODEL STATES
# Differential equations PK model
d/dt(A_centr)  = - CL/Vc*A_centr - Q/Vc*A_centr + A_periph*Q/Vp + INPUT1
d/dt(A_periph) =                   Q/Vc*A_centr - A_periph*Q/Vp

# Differential equations Friberg model
d/dt(A_prol)   = KTR*A_prol*(1-SLOPU*CONC)*(CIRC0/A_circ)^GAM - KTR*A_prol
d/dt(A_tr1)    = KTR*A_prol - KTR*A_tr1
d/dt(A_tr2)    = KTR*A_tr1 - KTR*A_tr2
d/dt(A_tr3)    = KTR*A_tr2 - KTR*A_tr3
d/dt(A_circ)   = KTR*A_tr3 - KTR*A_circ

# Initial conditions
A_prol(0)      = CIRC0
A_tr1(0)       = CIRC0
A_tr2(0)       = CIRC0
A_tr3(0)       = CIRC0
A_circ(0)      = CIRC0
********** MODEL PARAMETERS
# Define parameters
CIRC0          =  7.21   # Baseline circulating neutrophils (10^9 cells/L)
MTT            =  124    # Mean transit time (hours)
SLOPU          =  28.9   # Slope of drug effect (mL/ug)
GAM            = 0.239   # Exponent (.)

# Define regression parameters
Vp             = 1    # Peripheral volume (L)
Vc             = 1    # Central volume (L)
CL             = 1    # Clearance (L/hour)
Q              = 1    # Intercompartmental clearance (L/hour)
********** MODEL VARIABLES
# Calculate central concentration
CONC           = A_centr/Vc
# PD parameters
KTR            = 4/MTT
# Output for fitting
OUTPUT1        = A_circ  # (10^9 cells/L) Circulating neutrophils
********** MODEL REACTIONS
********** MODEL FUNCTIONS
********** MODEL EVENTS

### 9.6.3 Novak-Tyson Cell-Cycle

********** MODEL NAME
Novak-Tyson Model
********** MODEL NOTES
Novak-Tyson cell cycle model, described in J. theor. Biol. (1998) 195, 69-85
Used in IQRsim as an example model for quick model generation and testing.
********** MODEL STATES
% Comment 1
d/dt(Cyclin) = R1-R2-R3
d/dt(YT) = R4-R5-R6-R7+R8+R3
d/dt(PYT) = R5-R8-R9-R10+R11
d/dt(PYTP) = R12-R11-R13-R14+R9
d/dt(MPF) = R6-R4-R12-R15+R13
d/dt(Cdc25P) = R16
d/dt(Wee1P) = R17
d/dt(IEP) = R18
d/dt(APCstar) = R19

# Comment 2
Cyclin(0) = 0.0172
YT(0) = 0.0116
PYT(0) = 9e-04
PYTP(0) = 0.0198
MPF(0) = 0.073
Cdc25P(0) = 0.95
Wee1P(0) = 0.95
IEP(0) = 0.242
APCstar(0) = 0.3132
********** MODEL PARAMETERS
Ka = 0.1     % Comment 3
Kb = 1       # Comment 4
Kc = 0.01
Kd = 1
Ke = 0.1
Kf = 1
Kg = 0.01
Kh = 0.01
k1 = 0.01
k3 = 0.5
V2p = 0.005
V2pp = 0.25
V25p = 0.017
V25pp = 0.17
Vweep = 0.01
Vweepp = 1
kcak = 0.64
kpp = 0.004
kas = 2
kbs = 0.1
kcs = 0.13
kds = 0.13
kes = 2
kfs = 0.1
kgs = 2
khs = 0.15
********** MODEL VARIABLES
k2 = V2p+APCstar*(V2pp-V2p)
kwee = Vweepp+Wee1P*(Vweep-Vweepp)
k25 = V25p+Cdc25P*(V25pp-V25p)
********** MODEL REACTIONS
R1 = k1
R2 = k2*Cyclin
R3 = k3*Cyclin
R4 = kpp*MPF
R5 = kwee*YT
R6 = kcak*YT
R7 = k2*YT
R8 = k25*PYT
R9 = kcak*PYT
R10 = k2*PYT
R11 = kpp*PYTP
R12 = kwee*MPF
R13 = k25*PYTP
R14 = k2*PYTP
R15 = k2*MPF
R16 = kas*MPF*(1-Cdc25P)/(1+Ka-Cdc25P)-kbs*Cdc25P/(Kb+Cdc25P) {reversible}
R17 = kes*MPF*(1-Wee1P)/(1+Ke-Wee1P)-kfs*Wee1P/(Kf+Wee1P) {reversible}
R18 = kgs*MPF*(1-IEP)/(1+Kg-IEP)-khs*IEP/(Kh+IEP) {reversible}
R19 = kcs*IEP*(1-APCstar)/(1+Kc-APCstar)-kds*APCstar/(Kd+APCstar) {reversible}
********** MODEL FUNCTIONS
********** MODEL EVENTS

### 9.6.4 Parasitemia PK/PD

********** MODEL NAME
modelPKPD
********** MODEL NOTES
********** MODEL STATES
# PK model (dosing: INPUT1)
d/dt(Ac)    =  ka*Ad - CL/Vc*Ac - Q1/Vc*Ac + Q1/Vp1*Ap1
d/dt(Ap1)   =                     Q1/Vc*Ac - Q1/Vp1*Ap1

# PD model
d/dt(PD)    = GR - EMAX*Cc^hill/(EC50^hill+Cc^hill)
PD(0)       = PDbase
********** MODEL PARAMETERS
# PK parameters
Fabs1       = 1
ka          = 0.033
CL          = 0.192
Vc          = 0.474
Q1          = 1
Vp1         = 2

# PD parameters
PDbase      = 0
GR          = 0.05        # (1/hr)
EMAX        = 0.10        # (1/hr)
EC50        = 0.10        # (ug/ml)
hill        = 3           # (.)
********** MODEL VARIABLES
# Drug concentration (ug/mL)
Cc          = Ac/Vc
# Output
OUTPUT1     = Cc
OUTPUT2     = PD
********** MODEL REACTIONS
********** MODEL FUNCTIONS
********** MODEL EVENTS

### 9.6.5 Bouncing ball

********** MODEL NAME
Bouncing ball
********** MODEL NOTES
Simple example, implementing ODEs simulating a bouncing ball, using events.
********** MODEL STATES
d/dt(height) = velocity
d/dt(velocity) = -9.8

height(0) = 0
velocity(0) = 10
********** MODEL PARAMETERS
********** MODEL VARIABLES
********** MODEL REACTIONS
********** MODEL FUNCTIONS
********** MODEL EVENTS
Event1 = lt(height,0),height,0,velocity,-0.9*velocity

### 9.6.6 C-Functions

********** MODEL NAME
ModelC1
********** MODEL NOTES
IQRmodels can contain an additional section that allows to
include relatively arbitrary C code into models.
Some limitations of course exist. Generation of symbolic sensitivities might not
work, conversion to NONMEM and MONOLIX will not be possible.
********** MODEL STATES
d/dt(state1) = -k1*state1
d/dt(state2) =  k1*state1

state1(0) = 10
********** MODEL PARAMETERS
k1 = 0.5
********** MODEL VARIABLES
AA = test1()
BB = test2(state1,state2)
********** MODEL REACTIONS
********** MODEL FUNCTIONS
********** MODEL EVENTS
********** MODEL C FUNCTIONS
double test1() {
Rprintf("Hello World %f\n",2.3);
return(1);
}

double test2(double x, double y) {
if (x > y) return(x);
return(y);
}

### 9.6.7 Fantasy events

********** MODEL NAME
Fantasy Event Model
********** MODEL NOTES
Fantasy model based on the bouncing ball example but adding more events
testing the use of state contraints and changing both states and parameters
using events.
********** MODEL STATES
d/dt(height) = velocity
d/dt(velocity) = -9.8*p2
d/dt(X1) = 1

height(0) = height0
velocity(0) = 10
X1(0) = 0
********** MODEL PARAMETERS
p1 = 1
p2 = 1
height0 = 10
********** MODEL VARIABLES
v1 = mult(p1,p2)
********** MODEL REACTIONS
********** MODEL FUNCTIONS
mult(x,y) = x*y
********** MODEL EVENTS
Event1 = lt(height,-0.001),velocity,-0.9*velocity,height,0,p1,-p1
Event2 = ge(time,10),height,height+5,velocity,0,p1,3
Event3 = ge(time,15),p2,10

### 9.6.8 Novak-Tyson biochemical

********** MODEL NAME
Novak-Tyson Model
********** MODEL NOTES
Novak-Tyson cell cycle model, described in J. theor. Biol. (1998) 195, 69-85
********** MODEL STATE INFORMATION
Cyclin(0) = 0.0172
YT(0) = 0.0116
PYT(0) = 9e-04
PYTP(0) = 0.0198
MPF(0) = 0.073
Cdc25P(0) = 0.95
Wee1P(0) = 0.95
IEP(0) = 0.242
APCstar(0) = 0.3132
********** MODEL PARAMETERS
Ka = 0.1
Kb = 1
Kc = 0.01
Kd = 1
Ke = 0.1
Kf = 1
Kg = 0.01
Kh = 0.01
k1 = 0.01
k3 = 0.5
V2p = 0.005
V2pp = 0.25
V25p = 0.017
V25pp = 0.17
Vweep = 0.01
Vweepp = 1
kcak = 0.64
kpp = 0.004
kas = 2
kbs = 0.1
kcs = 0.13
kds = 0.13
kes = 2
kfs = 0.1
kgs = 2
khs = 0.15
********** MODEL VARIABLES
k2 = V2p+APCstar*(V2pp-V2p)
kwee = Vweepp+Wee1P*(Vweep-Vweepp)
k25 = V25p+Cdc25P*(V25pp-V25p)
********** MODEL REACTIONS
=> Cyclin : R1
vf = k1

Cyclin =>  : R2
vf = k2*Cyclin

Cyclin => YT : R3
vf = k3*Cyclin

MPF => YT : R4
vf = kpp*MPF

YT => PYT : R5
vf = kwee*YT

YT => MPF : R6
vf = kcak*YT

YT =>  : R7
vf = k2*YT

PYT => YT : R8
vf = k25*PYT

PYT => PYTP : R9
vf = kcak*PYT

PYT =>  : R10
vf = k2*PYT

PYTP => PYT : R11
vf = kpp*PYTP

MPF => PYTP : R12
vf = kwee*MPF

PYTP => MPF : R13
vf = k25*PYTP

PYTP =>  : R14
vf = k2*PYTP

MPF =>  : R15
vf = k2*MPF

<=> Cdc25P : R16
vf = kas*MPF*(1-Cdc25P)/(1+Ka-Cdc25P)
vr = kbs*Cdc25P/(Kb+Cdc25P)

<=> Wee1P : R17
vf = kes*MPF*(1-Wee1P)/(1+Ke-Wee1P)
vr = kfs*Wee1P/(Kf+Wee1P)

<=> IEP : R18
vf = kgs*MPF*(1-IEP)/(1+Kg-IEP)
vr = khs*IEP/(Kh+IEP)

<=> APCstar : R19
vf = kcs*IEP*(1-APCstar)/(1+Kc-APCstar)
vr = kds*APCstar/(Kd+APCstar)
********** MODEL FUNCTIONS
********** MODEL EVENTS


text:
Volume 316 - XXVI International Workshop on Deep-Inelastic Scattering and Related Subjects (DIS2018) - WG1: Structure Functions and Parton Densities
Strange and non-strange distributions from the collider data
S. Alekhin,* J. Bluemlein, S. Kulagin, S.O. Moch, R. Petti
*corresponding author
Full text: pdf
Pre-published on: 2018 September 20
Published on: 2018 November 23
Abstract
We check the stability of the ABMP16 fit with respect to
modifications of quark PDFs suggested in the recent literature: the strange sea enhancement and a positive non-vanishing $d/u$ ratio at $x \to 1$. These possibilities are examined using test versions of the ABMP16 PDF fit which demonstrate no need of those changes. Furthermore, we localize peculiar features in other analyses which are responsible for a different behaviour of the PDFs obtained.
The strange sea enhancement can be explained by a choice of the
PDF shapes being not flexible enough and therefore leading to an over-suppressed $d$-quark distribution. This suppression has to be compensated by a corresponding rise by the $s$-quark distribution.
As a result, an unusually large strange sea suppression factor
is obtained. The non-vanishing value of $d/u\vert _{x=1}$ becomes consistent with zero in case the higher-order QCD corrections to the D{\O} $W$-asymmetry data, which drive its value, are taken into account. Furthermore, the related $e$-asymmetry sample prefers a slightly negative value of $d/u\vert _{x=1}$, although it is consistent with zero. These clarifications support confidence in the PDF shapes used in the ABMP16 analyses.
DOI: https://doi.org/10.22323/1.316.0008
Open Access
Copyright owned by the author(s) under the term of the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License.