Mohammad Ehsanul Karim
2007-Apr-11 13:05 UTC
[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.
Kuhn, Max
2007-Apr-11 13:54 UTC
[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
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 at 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 messag...{{dropped}}