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

Wed Apr 11 15:05:35 CEST 2007

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
answer of the following questions:

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.

Thank you for your time.
Thanks in advance.

Mohammad Ehsanul Karim
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.