Danny - 3 months ago 12

R Question

In the package CoinMinD I ran the GM function on the following code:

`GM(c(2,4,5,1),0.05)`

The results are as follows:

`Original Intervals`

Lower Limit

[1] 0.03470190 0.10880765 0.15542939 0.01035233

Upper Limit

[1] 0.5266663 0.6718765 0.7349127 0.4413580

Adjusted Intervals

Lower Limit

[1] 0.03470190 0.10880765 0.15542939 0.01035233

Upper Limit

[1] 0.5266663 0.6718765 0.7349127 0.4413580

Volume

[1] 0.06918613

I am trying to create a matrix where the first two lines Lower Limit and the Upper Limit, appear as follows:

`[,1] [,2]`

[1,] 0.03470190 0.5266663

[2,] 0.10880765 0.6718765

[3,] 0.15542939 0.7349127

[4,] 0.01035233 0.4413580

Unforunately, when I save the code to an object:

`c<-GM(c(2,4,5,1),0.05)`

and run

`c[1]`

or

`GM(c(2,4,5,1),0.05) [1]`

I only get the last line of the results:

`> c[1]`

[1] 0.06918613

Answer

Interesting question! Luckily the function `GM`

is not a long one so we can look inside:

```
function (inpmat, alpha)
{
## ... those lines for computing results are left out...
cat("Original Intervals\n")
cat("Lower Limit\n")
print(GM.LL)
cat("Upper Limit\n")
print(GM.UL)
cat("Adjusted Intervals\n")
cat("Lower Limit\n")
print(LLA)
cat("Upper Limit\n")
print(ULA)
cat("Volume\n")
print(VOL)
}
```

All values you want to extract are only printed rather than exported. We can easily modify this function to export them.

```
myGM <- function (inpmat, alpha)
{
k = length(inpmat)
s = sum(inpmat)
chi = qchisq(1 - (alpha/k), df = 1)
pi = inpmat/s
GM.UL = (chi + 2 * inpmat + sqrt(chi * chi + 4 * inpmat *
chi * (1 - pi)))/(2 * (chi + s))
GM.LL = (chi + 2 * inpmat - sqrt(chi * chi + 4 * inpmat *
chi * (1 - pi)))/(2 * (chi + s))
LLA = 0
ULA = 0
for (r in 1:length(inpmat)) {
if (GM.LL[r] < 0)
LLA[r] = 0
else LLA[r] = GM.LL[r]
if (GM.UL[r] > 1)
ULA[r] = 1
else ULA[r] = GM.UL[r]
}
diA = ULA - LLA
VOL = round(prod(diA), 8)
## the lines for printing are removed
## replace it with the line below
list(limit = matrix(c(GM.LL, GM.UL, LLA, ULA), ncol = 4L,
dimnames = list(NULL, c("ori_low", "ori_up", "adj_low", "adj_up"))),
volume = VOL)
}
## test
myGM(c(2,4,5,1),0.05)
#$limit
# ori_low ori_up adj_low adj_up
#[1,] 0.03470190 0.5266663 0.03470190 0.5266663
#[2,] 0.10880765 0.6718765 0.10880765 0.6718765
#[3,] 0.15542939 0.7349127 0.15542939 0.7349127
#[4,] 0.01035233 0.4413580 0.01035233 0.4413580
#
#$volume
#[1] 0.06918613
```

You can easily extract the column you want from `$limit`

.