# [R] Programming Problem (for loop, random # control, 3 dimentional graph)

Kuhn, Max Max.Kuhn at pfizer.com
Wed Apr 11 16:14:29 CEST 2007

```Sorry, there was a typo in this line:

function(x) ks.test(x, "pnorm", mean = mean(x), sd = sd(x))\$p.value

Max

-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Kuhn, Max
Sent: Wednesday, April 11, 2007 9:55 AM
To: wildscop at yahoo.com; r-help at stat.math.ethz.ch
Subject: Re: [R] Programming Problem (for loop, random # control,3
dimentional graph)

Wow.

Some suggestions are below.

1. You should "vectorize" the for loop in the function g. I would do
something like this:

z1 <- rnorm(20 * 100, 0, 1)

Q1=0+1*z1*(1+.83*(1-exp(-G*z1))/(1+exp(-G*z1)))*(1+z1**2)**K

Q1Matrix <- matrix(Q1, ncol = 100)

ksStat <- function(x) ks.test(x, "pnorm", mean = mean(x), sd
=s(x))\$p.value

pValues <- apply(Q1Matrix, 2, ksStat)

p <- ifelse(pValues <= 0.05, 1, 0)

(please check the details. I have no idea of the purpose of this code).

Also, you might want to look at the outer function to compute across the
values of G and K if you want the results in a matrix (to avoid your
other loop). If you want the results in a data frame, you could use
expand.grid to create a grid of G and K combinations, then apply the
function to the rows (via apply).

2. I always put set.seed just prior the first function that call the
rng. Others probably have better advice.

3. To use traditional graphics like image our contour, save the data in
a matrix. For lattice graphics (my preference), the data frame approach
in #1 above is better. I think that scatterplot3d wants the results in
vector format (not matrix).

Good luck,

Max

-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Mohammad Ehsanul
Karim
Sent: Wednesday, April 11, 2007 9:06 AM
To: r-help at stat.math.ethz.ch
Subject: [R] Programming Problem (for loop, random # control,3
dimentional graph)

Dear List,

This is just a programming problem which i cannot seem
to figure out. I am trying to get a set of power from
a test (say, kolmogorov smirnov) out of a distribution
(say, G-K distribution) as follows. I am trying to
reduce to pain of writing the whole set of data points
(p# below) using "for" loop. However, I seem to have
some problem in it as the output "M" does not match
for the original and reduced form. I need to get

1. Is there any fault in the "for" loop in the Reduced
form?
2. To get exactly same result M in the next run, where
should i put set.seed()?
3. How to plot such plots in R having dimentions G, K
and vectorized M? In matlab, mesh(K,G,M) does the
trick. scatterplot3d(K,G,M) gives me error message.

wildscop at yahoo dot com
Institute of Statistical Research and Training
University of Dhaka

##########################################
#### Original form #######################

library(stats)
g=function(G,K,z1){
p=rep(0,100)
# 100 replication
for(i in 1:100){
z1=rnorm(20,0,1)
# standard normal with 20 sample size
Q1=0+1*z1*(1+.83*(1-exp(-G*z1))/(1+exp(-G*z1)))*(1+z1**2)**K
# density of G-K distribution
ks=ks.test(Q1, "pnorm", mean = mean(Q1), sd =
sqrt(var(Q1)))
pv=ks\$p.value
if(pv<=0.05) {p[i]=1} else {p[i]=0}
}
m.p=mean(p)
return(m.p)
}

# putting value of G -5,-4,-3,-2,-1,-.5,0,.5,1,2,3,4,5
# putting value of K -.5,0,.5,1,2,3,4,5
# in the function g()
p1=g(-5,-.5)
p2=g(-5,0)
p3=g(-5,0.5)
p4=g(-5,1)
p5=g(-5,2)
p6=g(-5,3)
p7=g(-5,4)
p8=g(-5,5)
p9=g(-4,-.5)
p10=g(-4,0)
p11=g(-4,0.5)
p12=g(-4,1)
p13=g(-4,2)
p14=g(-4,3)
p15=g(-4,4)
p16=g(-4,5)
p17=g(-3,-.5)
p18=g(-3,0)
p19=g(-3,.5)
p20=g(-3,1)
p21=g(-3,2)
p22=g(-3,3)
p23=g(-3,4)
p24=g(-3,5)
p25=g(-2,-.5)
p26=g(-2,0)
p27=g(-2,.5)
p28=g(-2,1)
p29=g(-2,2)
p30=g(-2,3)
p31=g(-2,4)
p32=g(-2,5)
p33=g(-1,-.5)
p34=g(-1,0)
p35=g(-1,.5)
p36=g(-1,1)
p37=g(-1,2)
p38=g(-1,3)
p39=g(-1,4)
p40=g(-1,5)
p41=g(-0.5,-0.5)
p42=g(-0.5,0)
p43=g(-0.5,0.5)
p44=g(-0.5,1)
p45=g(-0.5,2)
p46=g(-0.5,3)
p47=g(-0.5,4)
p48=g(-0.5,5)
p49=g(0,-0.5)
p50=g(0,0)
p51=g(0,.5)
p52=g(0,1)
p53=g(0,2)
p54=g(0,3)
p55=g(0,4)
p56=g(0,5)
p57=g(0.5,-.5)
p58=g(0.5,0)
p59=g(0.5,.5)
p60=g(0.5,1)
p61=g(0.5,2)
p62=g(0.5,3)
p63=g(0.5,4)
p64=g(0.5,5)
p65=g(1,-.5)
p66=g(1,0)
p67=g(1,.5)
p68=g(1,1)
p69=g(1,2)
p70=g(1,3)
p71=g(1,4)
p72=g(1,5)
p73=g(2,-.5)
p74=g(2,0)
p75=g(2,.5)
p76=g(2,1)
p77=g(2,2)
p78=g(2,3)
p79=g(2,4)
p80=g(2,5)
p81=g(3,-.5)
p82=g(3,0)
p83=g(3,.5)
p84=g(3,1)
p85=g(3,2)
p86=g(3,3)
p87=g(3,4)
p88=g(3,5)
p89=g(4,-.5)
p90=g(4,0)
p91=g(4,.5)
p92=g(4,1)
p93=g(4,2)
p94=g(4,3)
p95=g(4,4)
p96=g(4,5)
p97=g(5,-.5)
p98=g(5,0)
p99=g(5,0.5)
p100=g(5,1)
p101=g(5,2)
p102=g(5,3)
p103=g(5,4)
p104=g(5,5)
Mp<-c(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15,p16,p17,p18,p19
,p20,p21,p22,p23,p24,p25,p26,p27,p28,p29,p30,p31,p32,p33,p34,p35,p36,p37
,p38,p39,p40,p41,p42,p43,p44,p45,p46,p47,p48,p49,p50,p51,p52,p53,p54,p55
,p56,p57,p58,p59,p60,p61,p62,p63,p64,p65,p66,p67,p68,p69,p70,p71,p72,p73
,p74,p75,p76,p77,p78,p79,p80,p81,p82,p83,p84,p85,p86,p87,p88,p89,p90,p91
,p92,p93,p94,p95,p96,p97,p98,p99,p100,p101,p102,p103,p104)
M<-matrix(Mp,nrow=13,ncol=8,byrow=T)
M
##########################################
################ Result ##################
> M
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 0.33 0.23 0.47 0.60 0.85 0.97 0.98 0.97
[2,] 0.31 0.28 0.40 0.46 0.90 0.92 0.99 0.99
[3,] 0.28 0.26 0.34 0.61 0.85 0.95 0.99 0.99
[4,] 0.15 0.15 0.24 0.53 0.78 0.96 0.99 0.99
[5,] 0.04 0.04 0.19 0.35 0.84 0.94 0.96 0.99
[6,] 0.01 0.00 0.06 0.24 0.80 0.93 0.98 0.99
[7,] 0.00 0.00 0.02 0.27 0.73 0.94 0.97 0.98
[8,] 0.01 0.00 0.05 0.25 0.79 0.87 0.98 0.99
[9,] 0.07 0.03 0.22 0.44 0.71 0.94 0.97 1.00
[10,] 0.13 0.15 0.28 0.46 0.76 0.89 0.97 0.98
[11,] 0.27 0.27 0.35 0.55 0.86 0.97 0.99 1.00
[12,] 0.36 0.22 0.48 0.62 0.88 0.97 0.95 0.98
[13,] 0.41 0.29 0.34 0.59 0.88 0.96 0.99 0.99

##########################################
#### Reduced form ########################
library(stats)
g=function(G,K,r1){
p=rep(0,100)
for(i in 1:100){
z1=rnorm(r1,0,1)
Q1=0+1*z1*(1+.83*(1-exp(-G*z1))/(1+exp(-G*z1)))*(1+z1**2)**K
ks=ks.test(Q1, "pnorm", mean = mean(Q1), sd =
sqrt(var(Q1)))
pv=ks\$p.value
if(pv<=0.05) {p[i]=1} else {p[i]=0}
}
m.p=mean(p)
return(m.p)
}

G<-c(-5,-4,-3,-2,-1,-.5,0,.5,1,2,3,4,5)
K<-c(-.5,0,.5,1,2,3,4,5)
lg<-length(G)
lk<-length(K)
M<-matrix(rep(NA,lg*lk),nrow=lg,ncol=lk,byrow=T)
for(i in G[1]:G[lg]){
for(j in K[1]:K[lk]){
M[i,j]<-g(i,j,20)
}
}
M
##########################################
################ Result ##################
> M
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 0.62 0.76 0.97 0.99   NA   NA   NA   NA
[2,] 0.68 0.88 0.94 0.98   NA   NA   NA   NA
[3,] 0.65 0.89 0.95 0.98   NA   NA   NA   NA
[4,] 0.79 0.92 0.97 0.98   NA   NA   NA   NA
[5,] 0.70 0.85 0.95 1.00   NA   NA   NA   NA
[6,] 0.59 0.91 0.93 0.99   NA   NA   NA   NA
[7,] 0.59 0.91 0.93 0.99   NA   NA   NA   NA
[8,] 0.59 0.91 0.93 0.99   NA   NA   NA   NA
[9,] 0.59 0.91 0.93 0.99   NA   NA   NA   NA
[10,] 0.59 0.91 0.93 0.99   NA   NA   NA   NA
[11,] 0.59 0.91 0.93 0.99   NA   NA   NA   NA
[12,] 0.59 0.91 0.93 0.99   NA   NA   NA   NA
[13,] 0.59 0.91 0.93 0.99   NA   NA   NA   NA

________________________________________________________________________
____________
Looking for earth-friendly autos?
Browse Top Cars by "Green Rating" at Yahoo! Autos' Green Center.

______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help