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: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Kuhn, Max Sent: Wednesday, April 11, 2007 9:55 AM To: [EMAIL PROTECTED]; r-help@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: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Mohammad Ehsanul Karim Sent: Wednesday, April 11, 2007 9:06 AM To: r-help@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 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. ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. ---------------------------------------------------------------------- LEGAL NOTICE\ Unless expressly stated otherwise, this\ messa...{{dropped}} ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.