# Simulation metamodeling with constraints

**R snippets**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Last week I have posted about using simulation metamodeling to verify results of analytical solution of the model. After posting it I realized that the solution presented there can be improved by using knowledge of simulation model structure.

The story is about calculation of the variance of the random variable given by formula:

*q*(

*p*-2+

*q*/

*N*)

where *p* and *N* are model parameters and *q* is a random variable that has binomial distribution with *N* Bernoulli trials and probability of success 2-*p*. As I have written the properly specified formula for the metamodel has the following form:

**~ (**p

**+**I

**(**p

**^**2

**)**

**+**I

**(**p

**^**3

**)**

**+**I

**(**p

**^**4))

*** (**N

**+**I

**(**1

**/**N

**))**

where v is an estimate of the variance. Such a model was discussed in my last post (it is described in more detail there so I omit the repetition in this post).

*p*=1 and

*p*=2 the variance of the random variable is equal to 0. This is because in such case

*q*has no variance (it is equal to 1 and 0 respectively).

- the weight of cases where
*p*equals 1 or 2 might be increased (so we use Weighted Least Squares) - a restriction on values of prediction in points where
*p*equals 1 or 2 can be made to ensure that it is exactly 0 (so we use Constrained Least Squares)

**(**plyr

**)**

**(**limSolve

**)**

**<-**

**function**

**(**param

**)**

**{**

**<-**rbinom

**(**1000, param

**$**N, 2

**–**param

**$**p

**)**

**(**v

**=**var

**((**param

**$**p

**–**2

**+**q

**/**param

**$**N

**)**

*****q

**))**

**}**

**(**1

**)**

**<-**expand.grid

**(**p

**=**seq

**(**1, 2, len

**=**101

**)**, N

**=**1

**:**10

**)**

**<-**ddply

**(**data.set, .

**(**p,N

**)**, sim,.progress

**=**“text”

**)**

**<-**coef

**(**lm

**(**v

**~ (**p

**+**I

**(**p

**^**2

**)**

**+**I

**(**p

**^**3

**)**

**+**I

**(**p

**^**4

**))**

*** (**N

**+**I

**(**1

**/**N

**))**, data

**=**data.set

**))**

**<-**ifelse

**(**data.set

**$**v

**==**0, 100000, 1

**)**

**<-**coef

**(**lm

**(**v

**~ (**p

**+**I

**(**p

**^**2

**)**

**+**I

**(**p

**^**3

**)**

**+**I

**(**p

**^**4

**))**

*** (**N

**+**I

**(**1

**/**N

**))**, data

**=**data.set,

**=**weights

**))**

**<-**with

**(**data.set, cbind

**(**1, p, p

**^**2, p

**^**3, p

**^**4, N, 1

**/**N,

*****N, p

**/**N, p

**^**2

*****N, p

**^**2

**/**N,

**^**3

*****N, p

**^**3

**/**N, p

**^**4

*****N, p

**^**4

**/**N

**))**

**(**A

**)**

**<-**names

**(**lm.coef

**)**

**<-**data.set

**$**v

**<-**A

**[**B

**==**0,

**]**

**<-**B

**[**B

**==**0

**]**

**<-**lsei

**(**A, B, E, F

**)$**X

**<-**c

**(**32,

**–**88, 88,

**–**38, 6,

**–**8,

**–**26, 20,

**–**18,

**–**79, 7, 36,

**–**1,

**–**6

**)**

**<-**cbind

**(**lm.coef, wlm.coef, lsei.coef

**)**

**–**true.coef

**(**coef.dev, digits

**=**2

**)**

**(**coef.dev

**[**,1

**]**, pch

**=**19

**)**

**(**coef.dev

**[**,2

**]**, 1

**:**nrow

**(**coef.dev

**)**, col

**=**“red”, pch

**=**19

**)**

**(**v

**=**0, col

**=**“gray”, lty

**=**2

**)**

- the precision of estimation is improved by using additional information on problem structure;
- the estimates from WLS and CLS are identical to the third decimal place.

**leave a comment**for the author, please follow the link and comment on their blog:

**R snippets**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.