hey, my function doesn?t work. can somebody help me?
the graphic doesn?t work and also the function. thnx a lot.
N=10
n=100
p_0=c(1/5,1-1/5)
power = function(p,m) {
set.seed(1000)
H=matrix(0,nrow=N,ncol=1)
for(i in 1:N) {
x <- matrix(rnorm(n, 0, 0.5), ncol = m)
y <- matrix(rnorm(n, 0, 0.8), ncol = m)
l <- diag(cor(x, y))
q_1 = qnorm(0.05, 0, 0.05)
q_2 = qnorm(1 - 0.05, 0, 0.05)
p <- (l^2)/sum(l^2)
H[i] <- sum(p_0*log(p_0)) - sum(p * log(p))
}
1- mean(q_1 <= H & H <= q_2)
}
m=seq(10,50,len=10)
f=outer(p,m,Vectorize(power))
persp(p,m,power,theta=-50,phi=30,d=4,border="black")
--
View this message in context:
http://r.789695.n4.nabble.com/the-function-doesn-t-work-tp2714105p2714105.html
Sent from the R help mailing list archive at Nabble.com.
The message is clear. Just resove this problem before posting a
terribly general and so not useful "it does not work".
Best
mario
> f=outer(p,m,Vectorize(power))
Error in outer(p, m, Vectorize(power)) : object 'p' not found
> persp(p,m,power,theta=-50,phi=30,d=4,border="black")
Error in persp(p, m, power, theta = -50, phi = 30, d = 4, border =
"black") :
object 'p' not found
On 26-Sep-10 08:39, jethi wrote:> hey, my function doesn?t work. can somebody help me?
> the graphic doesn?t work and also the function. thnx a lot.
>
> N=10
> n=100
>
> p_0=c(1/5,1-1/5)
>
> power = function(p,m) {
> set.seed(1000)
> H=matrix(0,nrow=N,ncol=1)
> for(i in 1:N) {
> x<- matrix(rnorm(n, 0, 0.5), ncol = m)
> y<- matrix(rnorm(n, 0, 0.8), ncol = m)
> l<- diag(cor(x, y))
>
> q_1 = qnorm(0.05, 0, 0.05)
> q_2 = qnorm(1 - 0.05, 0, 0.05)
> p<- (l^2)/sum(l^2)
> H[i]<- sum(p_0*log(p_0)) - sum(p * log(p))
> }
> 1- mean(q_1<= H& H<= q_2)
> }
> m=seq(10,50,len=10)
> f=outer(p,m,Vectorize(power))
> persp(p,m,power,theta=-50,phi=30,d=4,border="black")
--
Ing. Mario Valle
Data Analysis and Visualization Group | http://www.cscs.ch/~mvalle
Swiss National Supercomputing Centre (CSCS) | Tel: +41 (91) 610.82.60
v. Cantonale Galleria 2, 6928 Manno, Switzerland | Fax: +41 (91) 610.82.82
hi, sorry but i can?t remove the problem.but i change the programm a little
bit. i didn?t work with r programm before, so its really hard for me to find
my problems. :)
N=5
n=100
p_0=c(1/5,1-1/5)
power = function(k1) {
set.seed(1000)
H=matrix(0,nrow=N,ncol=1)
for(i in 1:N) {
x <- matrix(rnorm(n, 0, 0.5), ncol =m1)
y <- matrix(rnorm(n, 0, 0.8), ncol = m1)
l <- diag(cor(x, y))
q_1 = qnorm(0.05, 0, 0.05)
q_2 = qnorm(1 - 0.05, 0, 0.05)
p <- (l^2)/sum(l^2)
H[i] <- sum(p_0*log(p_0)) - sum(p * log(p))
}
1- mean(q_1 <= H & H <= q_2)
}
m1=seq(0,n/2,len=10)
k1=1/m1
output <- power(k1)
f=outer(k1,Vectorize(power))
--
View this message in context:
http://r.789695.n4.nabble.com/the-function-doesn-t-work-tp2714105p2714123.html
Sent from the R help mailing list archive at Nabble.com.
This is worse than before and getting pretty silly. Now you are calling outer
with a function that only takes one argument.
R might be hard for you, but mind reading is even harder for most of us. To get
help you need to explain clearly and sensibly what it is you want to do.
Look at your code. You generate (the same) random numbers every time and do not
use them. diag(cor(x,y)) is always a vector of ones. So p is the same at every
stage of the loop as well. You put q_1 and q_2 inside the loop, but only use
the outside. This is not R, this is just a failure to think clearly about what
you are doing.
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On Behalf Of jethi
Sent: Sunday, 26 September 2010 5:29 PM
To: r-help at r-project.org
Subject: Re: [R] the function doesn?t work
hi, sorry but i can?t remove the problem.but i change the programm a little
bit. i didn?t work with r programm before, so its really hard for me to find
my problems. :)
N=5
n=100
p_0=c(1/5,1-1/5)
power = function(k1) {
set.seed(1000)
H=matrix(0,nrow=N,ncol=1)
for(i in 1:N) {
x <- matrix(rnorm(n, 0, 0.5), ncol =m1)
y <- matrix(rnorm(n, 0, 0.8), ncol = m1)
l <- diag(cor(x, y))
q_1 = qnorm(0.05, 0, 0.05)
q_2 = qnorm(1 - 0.05, 0, 0.05)
p <- (l^2)/sum(l^2)
H[i] <- sum(p_0*log(p_0)) - sum(p * log(p))
}
1- mean(q_1 <= H & H <= q_2)
}
m1=seq(0,n/2,len=10)
k1=1/m1
output <- power(k1)
f=outer(k1,Vectorize(power))
--
View this message in context:
http://r.789695.n4.nabble.com/the-function-doesn-t-work-tp2714105p2714123.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
R-help at r-project.org 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.
It is not a problem of not knowing R.
It is a problem of reasoning. if you use m1 and not assign to it a value
beforehand it is difficult your function works.
And this will happen in any language, not only R.
Maybe explaining what you are trying to do helps. To do this try to add
comments (starting with #) to your code.
Reading error messages is always useful too.
Best
mario
On 26-Sep-10 09:28, jethi wrote:> hi, sorry but i can?t remove the problem.but i change the programm a little
> bit. i didn?t work with r programm before, so its really hard for me to
find
> my problems. :)
>
>
> N=5
> n=100
>
>
> p_0=c(1/5,1-1/5)
>
> power = function(k1) {
> set.seed(1000)
> H=matrix(0,nrow=N,ncol=1)
>
> for(i in 1:N) {
>
> x<- matrix(rnorm(n, 0, 0.5), ncol =m1)
> y<- matrix(rnorm(n, 0, 0.8), ncol = m1)
> l<- diag(cor(x, y))
>
>
>
> q_1 = qnorm(0.05, 0, 0.05)
> q_2 = qnorm(1 - 0.05, 0, 0.05)
> p<- (l^2)/sum(l^2)
>
> H[i]<- sum(p_0*log(p_0)) - sum(p * log(p))
>
> }
> 1- mean(q_1<= H& H<= q_2)
>
> }
> m1=seq(0,n/2,len=10)
> k1=1/m1
> output<- power(k1)
> f=outer(k1,Vectorize(power))
--
Ing. Mario Valle
Data Analysis and Visualization Group | http://www.cscs.ch/~mvalle
Swiss National Supercomputing Centre (CSCS) | Tel: +41 (91) 610.82.60
v. Cantonale Galleria 2, 6928 Manno, Switzerland | Fax: +41 (91) 610.82.82
i?m really sorry, once again. ok i will try to explain what i have to
programm. i want to programm a powerfunction.
i have to research if the correlations in a bivariate random sample are
homogeneous. for that i saperate the random sample in m blocks and
calculate the correlation of each block(partial sample). than i look at the
random variable p=cor(x1,y1)^2/sum(cor(x,y)^2) which is a probability. my
null hypothesis is that all correlations are homogeneous. if we have this
situation the entropy take the value log(m).
my test based on the entropy. so my teststatic is H=log(m)-sum(p*log(p)).
the following programm was actually, which i have worked the most of the
time. i hope that i don?t confuse u too much.
i of course i hope u understand my problem and my theme.
n=1000
m=2
k=n/m
N=100
myfun <- function(n, m, alpha = .05, seeder = 1000) {
l=matrix(0,nrow=m,ncol=N)
for(i in 1:N){
set.seed(i)
for(j in 1:m){
x=rnorm(n,0,0.5)
y=rnorm(n,0,0.8)
l[j,i]=cor((x[(((j-1)*k)+1):(((j-1)*k)+k)]),
(y[(((j-1)*k)+1):(((j-1)*k)+k)]))
}
}
gute <- function(x,m) {
q_1 <- qnorm(alpha, 0, 0.05)
q_2 <- qnorm(1 - alpha, 0, 0.05)
p=matrix(0,nrow=m,ncol=N)
H=matrix(0,nrow=N,ncol=1)
for(i in 1:N){
for (j in 1:m){
p[j,i]=x[j]^2/sum(x[,i]^2)
}
H[i]=log(m)-sum(p[,i]*log(p[,i]))
}
1 - mean(q_1 <= H & H <= q_2)
}
output <- gute(x = l,m=m)
return(output)
}
--
View this message in context:
http://r.789695.n4.nabble.com/the-function-doesn-t-work-tp2714105p2714137.html
Sent from the R help mailing list archive at Nabble.com.
Hi Jethi,
Please look at this code, and the graph that is created in the code.
This all leads me to suspect you are trying to use the wrong formula
at the end. Since there was some confusion with it, I just removed
the stuff with apply().
Josh
##################################
#### I moved stuff around to assign all global variables at the beginning
## Assign global variables; this works
N = 10
n = 100
m = 2
k = n/m
l = matrix(0, nrow=m, ncol=N)
p = matrix(0, nrow=m, ncol=N)
H = matrix(0, nrow=N, ncol=1)
alpha = 0.05
q_1 <- qnorm(alpha, 0, 0.05)
q_2 <- qnorm(1 - alpha, 0, 0.05)
## Calculate correlations; this works
for(i in 1:N){
x = rnorm(n,0,0.5)
y = rnorm(n,0,0.8)
for(j in 1:m){
l[j,i] = cor(
(x[(((j-1)*k)+1):(((j-1)*k)+k)]),
(y[(((j-1)*k)+1):(((j-1)*k)+k)]))
}
}
## Calculate probabilities; this works
for(i in 1:N){
for (j in 1:m){
p[j,i] = (l[j,i]^2) / sum(l[,i]^2)
}
}
## Calculate this other thing; this works, but is not what you want
for(i in 1:N){
H[i]=log(m)-sum(p[,i]*log(p[,i]))
}
1 - mean(q_1 <= H & H <= q_2)
## lets break this down then, and see what each part does
## first we are going to redefine H. H = log(m) - zz
## where zz = sum(p[,i] * log(p[,i]))
## create the matrix zz
zz <- matrix(0, nrow=N, ncol=1)
for(i in 1:N) {
zz[i] <- sum(p[,i] * log(p[,i]))
}
## You have a logical test (q_1 <= H <= q_2)
## that is true when H is between q_1 and q_2
## Lets plot our data to see what is going on
plot(x = seq_along(zz), y = zz, ylim = c(-1, 2), col = "black", pch =
16)
points(x = seq_along(zz), y = H, col = "blue", pch = 16)
abline(h = q_1, col = "red")
abline(h = q_2, col = "red")
abline(h = log(m), col = "green")
## The black dots are the values of zz (sum(p[,i] * log(p[,i])))
## the blue dots are your actual values of H
## the red lines are the region selected by q_1 and q_2 (based on alpha)
## for your logical test to evaluate TRUE, the blue dots MUST be
within the red lines
## otherwise it will be FALSE
## finally, the green line is log(m), you can see that log(m) is so far outside
## of the selection region (between red lines), that it pushes
## every sum(p[,i] * log(p[,i])) outside
## To change this, you need to lower log(m) (~ m = .8 would work)
## or you need to increase sum(p[,i] * log(p[,i])), which is tricky
## because log(1) = 0, and the log of anything in the interval (0, 1)
## is a negative number. If p is a probability, it must be between [0, 1]
## which means you are guaranteed to be multiplying a positive number by
## a negative number, or 0. This in turn means the sum of all of these
## will itself be negative, or 0. So, since H = log(m) - zz, and we have
## established that zz must be 0 or negative, the minimum bound of H is log(m)
## and you can see from the graph (log(m) is the green line), that log(m)
## is far above your selection region
#############################
On Sun, Sep 26, 2010 at 5:13 PM, jethi <kartija at hotmail.com>
wrote:>
> hi josh, and really thnx again.
>
> i have now 2 problems . the first one ist if i take ur idea and programm
> like u:
>
> N = 10
> n = 100
> m = 2
> k = n/m
> l = matrix(0,nrow=m,ncol=N)
> p=matrix(0,nrow=m,ncol=N)
>
> alpha = 0.05
> q_1 <- qnorm(alpha, 0, 0.05)
> q_2 <- qnorm(1 -alpha, 0, 0.05)
> for(i in 1:N){
> ?x=rnorm(n,0,0.5)
> ?y=rnorm(n,0,0.8)
> ?for(j in 1:m){
> ? ?l[j,i] = cor(
> ? ? ? (x[(((j-1)*k)+1):(((j-1)*k)+k)]),
> ? ? ? (y[(((j-1)*k)+1):(((j-1)*k)+k)]))
>
>
> ?}
> }
>
>
> for(i in 1:N){
> for (j in 1:m){
> p[j,i]=l[j,i]^2/sum(l[,i]^2)
> }
> }
>
> p2 <- apply(l, 2, function(z) {(z^2)/sum(z^2)})
> identical(p, p2)
>
> H=matrix(0,nrow=N,ncol=1)
> for(i in 1:N){
> H[i]=log(m)-sum(p[,i]*log(p[,i]))
> }
> H1<-apply(p2,2,function(t){(log(m)-sum(t*log(t)))})
> identical(H1, H)
>
>
> it gives me, that identical(H1, H) is False. is it because of H? H is a
> matrix and H1 is a list, isn?t it?
> ok if i leave H1 out and programm further to get my powerfunction, it gives
> me always the value 1. even i change my alpha or the number of n, m or N.:
>
>
> N = 10
> n = 100
> m = 2
> k = n/m
> l = matrix(0,nrow=m,ncol=N)
> p=matrix(0,nrow=m,ncol=N)
>
> alpha = 0.05
> q_1 <- qnorm(alpha, 0, 0.05)
> q_2 <- qnorm(1 -alpha, 0, 0.05)
> for(i in 1:N){
> ?x=rnorm(n,0,0.5)
> ?y=rnorm(n,0,0.8)
> ?for(j in 1:m){
> ? ?l[j,i] = cor(
> ? ? ? (x[(((j-1)*k)+1):(((j-1)*k)+k)]),
> ? ? ? (y[(((j-1)*k)+1):(((j-1)*k)+k)]))
>
>
> ?}
> }
> for(i in 1:N){
> for (j in 1:m){
> p[j,i]=l[j,i]^2/sum(l[,i]^2)
> }
> }
>
> p2 <- apply(l, 2, function(z) {(z^2)/sum(z^2)})
> identical(p, p2)
>
> H=matrix(0,nrow=N,ncol=1)
> for(i in 1:N){
> H[i]=log(m)-sum(p[,i]*log(p[,i]))
> }
> 1-mean(q_1<=H & H <=q_2)
>
> regards
> jethi
>
>
> --
> View this message in context:
http://r.789695.n4.nabble.com/the-function-doesn-t-work-tp2714105p2714815.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at r-project.org 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.
>
--
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.com/
wow thnx a lot josh. so now i understand the most of the things. so my power function is depent on the number of the blocks "m". how u said, if i take m=0.8 or lower i get a different value than a 1. ofcourse a number between 0 or 1. but i can just take positive interger number for m, because m representet die value for the blocks and they could be only a positiv number. moreover the graphic show that my H will never reach log(m) because like u said p are probilities. so if i would now plot the powerfunction for different m(positiv integer), the graph would be a constant value. am i right? thank u very very much. its help me a lot. -- View this message in context: http://r.789695.n4.nabble.com/the-function-doesn-t-work-tp2714105p2714875.html Sent from the R help mailing list archive at Nabble.com.
As long as p <= 1, the minimum value of H will be log(m). Here are
some more (I think clearer) graphs. They show the basic function
p[,i] * log(p[,i] (the log function defaults to natural logarithm),
for values ranging from 0 to 1. Then include log(m) for different
values of m. I included the code, and also attached a PDF in incase
there was any trouble running the code.
x <- seq(0, 1, by = .01)
par(mfrow = c(3, 2))
plot(x = x, y = log(x ^ x), type = "l",
ylab = bquote(f(x) == ln(x^x)), main = "Basic graph")
plot(x = x, y = -log(x ^ x), type = "l",
ylab = bquote(f(x) == -ln(x^x)), main = "Basic negative graph")
plot(x = x, y = log(0.5) - log(x ^ x), type = "l",
ylab = bquote(f(x) == ln(0.5) - ln(x^x)), main = "m = 0.5")
plot(x = x, y = log(1) - log(x ^ x), type = "l",
ylab = bquote(f(x) == ln(1) - ln(x^x)), main = "m = 1")
plot(x = x, y = log(2) - log(x ^ x), type = "l",
ylab = bquote(f(x) == ln(2) - ln(x^x)), main = "m = 2")
plot(x = x, y = log(6) - log(x ^ x), type = "l",
ylab = bquote(f(x) == ln(6) - ln(x^x)), main = "m = 6")
So the way I see it, you have 3 ways to get different values:
1) Increase your matrix p
2) Decrease m
3) Increase your selection region (this in turn depends on alpha, the
mean, and the standard deviation)
Cheers,
Josh
On Sun, Sep 26, 2010 at 7:58 PM, jethi <kartija at hotmail.com>
wrote:>
> wow thnx a lot josh. so now i understand the most of the things. so my
power
> function is depent on the number of the blocks "m". how u said,
if i take
> m=0.8 or lower i get a different value than a 1. ofcourse ?a number between
> 0 or 1. but i can just take positive interger number for m, because m
> representet die value for the blocks and they could be only a positiv
> number. moreover the graphic show that my H will never reach log(m) because
> like u said p are probilities. so if i would now plot the powerfunction for
> different m(positiv integer), the graph would be a constant value. am i
> right?
>
> thank u very very much. its help me a lot.
> --
> View this message in context:
http://r.789695.n4.nabble.com/the-function-doesn-t-work-tp2714105p2714875.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at r-project.org 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.
>
--
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.com/
-------------- next part --------------
A non-text attachment was scrubbed...
Name: graphs.pdf
Type: application/pdf
Size: 24875 bytes
Desc: not available
URL:
<https://stat.ethz.ch/pipermail/r-help/attachments/20100926/c8ecce9d/attachment.pdf>