Commit 6f7e253c by Eric Coissac

Delete Bricolage.Rmd

parent d955695d
---
title: "R Notebook"
output:
html_document:
df_print: paged
---
```{r}
Xs = test2
nperm=100
```
```{r}
ix = c(4,3,6,5,8,7)
iy = 1
r = corls(Xs,permutations = how(nperm = nperm))
rx = r[ix,ix]
ry = r[iy,ix]
rxi = r[ix[(ix %% 2)>0],ix[(ix %% 2)>0]]
ryi = r[iy,ix[(ix %% 2)>0]]
```
```{r}
r2 = (solve(rx) %*% ry) * ry
r2
```
```{r}
sum(r2)
```
```{r}
r2 = (solve(rxi) %*% ryi) * ryi
r2
```
```{r}
sum(r2)
```
```{r}
v = varls(Xs,permutations = how(nperm = nperm))
vt = v[c(iy,ix[(ix %% 2)>0]),c(iy,ix[(ix %% 2)>0])] * (nrow(Xs)-1)
vr = v[c(iy+1,ix[(ix %% 2)==0]),c(iy,ix[(ix %% 2)>0])] * (nrow(Xs)-1)
vt
vr
```
```{r}
#diag(vr)=0
vp = vt - vr
vp
```
```{r}
d = sqrt(diag(vp))
n = d %o% d
r = vp / n
r
```
```{r}
rx = r[ix[(ix %% 2)==0]/2,ix[(ix %% 2)==0]/2]
ry = r[iy,ix[(ix %% 2)==0]/2]
r2 = (solve(rx) %*% ry) * ry
r2
```
```{r}
sum(r2)
```
```{r}
rp = solve(r)
rp = rp / sqrt(diag(rp) %o% diag(rp))
rp
rp^2
```
```{r}
library(ProcMod)
<<<<<<< HEAD
r2s = seq(from=0.0, to = 0.9, by = 0.1)
n.sim=1000
r2.obs = matrix(0,nrow = n.sim,ncol = length(r2s))
d1 = simulate_matrix(10,2)
for (j in seq_along(r2s)) {
#print(c("r2",as.character(r2s[j])))
for (i in seq_len(n.sim)) {
# print(c("round",as.character(i)))
d2 = simulate_correlation(d1,1,r2s[j])
#r2.obs[i,j] = cor5(d1,d2)^2
r2.obs[i,j] = ((max(0,.covls(d1,d2)-.rcovls(d1,d2)))/(sqrt(.covls(d1,d1)-.rcovls(d1,d1))*sqrt(.covls(d2,d2)-.rcovls(d2,d2))))^2
# v = varls(d1,d2,permutations = 10)
# s = sqrt(diag(v))
# r2.obs[i,j] = (v / outer(s,s))[1,2]^2
=======
r2s = seq(from=0.1, to = 0.9, by = 0.1)
n.sim=20
r2.obs = matrix(0,nrow = n.sim,ncol = length(r2s))
d1 = simulate_matrix(10,2000)
for (j in seq_along(r2s)) {
print(c("r2",as.character(r2s[j])))
for (i in seq_len(n.sim)) {
print(c("round",as.character(i)))
d2 = simulate_correlation(10,3000,d1,r2s[j])
v = varls2(d1,d2,permutations = 10)
s = sqrt(diag(v))
r2.obs[i,j] = (v / outer(s,s))[1,2]^2
>>>>>>> ccb4d408c4936c553b8d097d94022285a7bc89a7
}
}
boxplot(r2.obs)
cbind(R2_theo=r2s,
R2_obs=colMeans(r2.obs),
R2_sd=apply(r2.obs,MARGIN = 2,FUN = sd))
```
<<<<<<< HEAD
```{r}
a= simulate_matrix(40,1)
b= simulate_matrix(40,2)
x = varls(a,b,permutations = NULL)
y = corls(a,b,permutations = 100)
y^2
m = x - y
#(sqrt(outer(diag(x),diag(x)))/m)^2
#sqrtm(outer(diag(x),diag(x)))
```
5 -> 6
10 -> 12 6
20 -> 28 14
30 -> 44 22
40 -> 60 30 20
```{r}
i=c(5,10,20,30,40)
j=c(6,12,28,44,60)
lm(j~i)
plot(i,j,xlim=c(0,50),ylim = c(-10,70))
abline(lm(j~i))
plot(sqrt(i),sqrt(j))
lm(sqrt(j)~sqrt(i))
plot(i^2,j^2)
lm(i^2~j^2)
```
```{r}
a = simulate_matrix(20,1)
b = vector(mode = "list",9)
for (i in seq_along(b)) {
b[[i]] = simulate_correlation(a,1,i/10)
}
b = c(list(a),b)
var(do.call(cbind,b))
var2(do.call(cbind,b))
v.perm=varls(as.procmod.frame(b),permutations = 100,rcovls = TRUE)
v.perm
m = apply(attributes(v.perm)$rcov,MARGIN = c(1,2),mean,na.rm=TRUE)^2
mean(m[upper.tri(m,diag = TRUE)],na.rm=TRUE)
```
```{r}
library(ProcMod)
ps = c(seq(1,10,by=1),seq(15,100,by=5),1000)
qs = ps
ns = seq(5,20, by=5)
n.sim=1000
r2.obs = array(0,dim = c(length(ns),length(ps),length(qs),n.sim))
for (ni in seq_along(ns)) {
N=ns[ni]
for (pi in seq_along(ps)) {
P = ps[pi]
for (qi in seq_along(qs)) {
Q = qs[qi]
for (si in 1:n.sim) {
not_ok = TRUE
while(not_ok) {
d1 = simulate_matrix(N,P,equal.var = TRUE)
d2 = simulate_matrix(N,Q,equal.var = TRUE)
r2=tryCatch({v = varls(d1,d2,permutations = 0)
s = sqrt(diag(v))
(v / outer(s,s))[1,2]^2
},
error=function(cond)
{#message("Redo simulation")
NA}
)
r2.obs[ni,pi,qi,si] = r2
not_ok=is.na(r2)
}
}
}
}
}
```
```{r}
R2_obs=apply(r2.obs,MARGIN = c(1,2,3),FUN = mean)
R2_sd=apply(r2.obs,MARGIN = c(1,2,3),FUN = sd)
dimnames(R2_obs) = list(ns,ps,qs)
N = rep(ns,length(ps)*length(qs))
P = rep(rep(ps,rep(length(ns),length(ps))),length(qs))
Q = rep(qs,rep(length(ns)*length(ps),length(qs)))
R2_obs.tab = tibble(N=N,
P=P,
Q=Q,
Pmin = mapply(min,P,Q),
Qmax = mapply(max,P,Q),
mean=as.numeric(R2_obs),
sd=as.numeric(R2_sd))
R2_obs.tab$NPmin = mapply(min,R2_obs.tab$Pmin,N)
R2_obs.tab$NPmax = mapply(max,R2_obs.tab$Pmin,N)
write_csv2(R2_obs.tab,"R2_H0.csv")
#p = rep(ps,rep(length(qs),length(ps)))
#q = rep(qs,length(ps))
#R2_obs.tn = data.frame(p,q,N=N,mean=as.numeric(R2_obs),sd=as.numeric(R2_sd))
```
```{r}
library(expm)
```
# Modle 1 : l'axe suivant se prend le reste
```{r}
R2.alea.predict.single = function(N,nP,nQ) {
# N : nombre de points
# nP : Nb de var dans petit tableau
# nQ : Nb de var dans grand tableau
# Axes principaux : P = (N-2)/(N-1) ; Q = 1/(N-1)
Axes.main = min(nP,(N-1))
i = 0:(Axes.main-1)
Q = 1/(N-1)
P = (N-2)/(N-1)
Pis = P**i
R2.main = sum(Pis)*Q
# Axes secondaires : Q = Axes.main/(N-1)^2 ; P = 1 - Q ;
Axes.comp = max(nP,nQ) - Axes.main
if (Axes.comp >0) {
i = 0:(Axes.comp-1)
Q = 1/(N-1)^2
P = 1 - Q
Pis = (P**i)
R2.main = P**Axes.comp*R2.main+Q*sum(Pis)
}
return(R2.main)
}
```
# Modle 2 : l'axe suivant se prend le reste
```{r}
R2.alea.predict.single = function(N,nP,nQ) {
# N : nombre de points
# nP : Nb de var dans petit tableau
# nQ : Nb de var dans grand tableau
# Axes principaux : P = (N-2)/(N-1) ; Q = 1/(N-1)
Axes.main = min(nP,(N-1))
Rres = numeric(Axes.main)
R2.main=0
for (i in 1:Axes.main) {
Rres[i]=1-R2.main
if (i==1) k = 0
else k = sum(Rres[1:(i-1)]^2)
# print(k / (N-1)^2 * Rres[i])
R2.main = R2.main + Rres[i] * ( 1/(N-1) - k/(N-1)^2)
}
# R2.main = R2.main - 2*Axes.main / (N-1)^2
# print(2*Axes.main / (N-1)^2)
# Axes secondaires : Q = Axes.main/(N-1)^2 ; P = 1 - Q ;
k=sum(Rres^2) + (1 - R2.main)^2
Axes.comp = max(nP,nQ) - Axes.main
if (Axes.comp >0) {
for (i in 1:Axes.comp) {
#print(R2.main)
R2.main = R2.main + (1-R2.main) * k/(N-1)^2
}
}
return(R2.main)
}
R2.alea.predict =Vectorize(R2.alea.predict.single,
vectorize.args = c("N","nP","nQ"))
```
# Modle 3 : l'axe suivant se prend le reste
```{r}
R2.alea.predict.single = function(N,nP,nQ) {
# N : nombre de points
# nP : Nb de var dans petit tableau
# nQ : Nb de var dans grand tableau
# Axes principaux : P = (N-2)/(N-1) ; Q = 1/(N-1)
Axes.main = min(nP,(N-1))
Rres = numeric(Axes.main)
R2.main=0
for (i in 1:Axes.main) {
Rres[i]=1-R2.main
R2.main = R2.main + Rres[i] * ( 1/(N-1) - (Axes.main - i )/(N-1)^2)
}
# R2.main = R2.main - 2*Axes.main / (N-1)^2
# print(2*Axes.main / (N-1)^2)
# Axes secondaires : Q = Axes.main/(N-1)^2 ; P = 1 - Q ;
k=sum(Rres) + (1 - R2.main)
Axes.comp = max(nP,nQ) - Axes.main
if (Axes.comp >0) {
for (i in 1:Axes.comp) {
#print(R2.main)
R2.main = R2.main + (1-R2.main) * k/(N-1)^2
}
}
return(R2.main)
}
R2.alea.predict =Vectorize(R2.alea.predict.single,
vectorize.args = c("N","nP","nQ"))
```
```{r}
library(tidyverse)
R2_obs.tab = read_csv2("R2_H0.csv")
```
```{r}
R2_obs.tab.PQ = R2_obs.tab %>% filter(Qmax==Pmin)
```
```{r}
with(data = R2_obs.tab.PQ,
plot(mean,
R2.alea.predict(N,Pmin,Qmax),
col=rainbow(4)[N/5],
pch=(1+(Qmax==Pmin)),
cex=0.5,
log="")
)
abline(b=1,a=0)
abline(h=1,lty=2)
abline(v=1,lty=2)
abline(lm(R2.alea.predict(N,Pmin,Qmax)~mean,data = R2_obs.tab.PQ))
#summary(lm(mean~R2.alea.predict(N,Pmin,Qmax),data = R2_obs.tab.PQ))
legend("bottomright",legend = (1:4)*5,fill=rainbow(4),cex = 0.5)
```
```{r}
with(data = R2_obs.tab.PQ,
plot(mean,
R2.alea.predict(N,Pmin,Qmax)-mean,
col=rainbow(4)[N/5],
pch=(1+(Qmax==Pmin)),
cex=0.5,
log="")
)
abline(h=1)
abline(v=2)
abline(v=c(5,10,15,20)*2.5,lty=2)
#abline(lm(log10(R2.alea.predict(N,Pmin,Qmax)/mean)[-1]~log10(Pmin)[-1],data=R2_obs.tab.PQ))
legend("topright",legend = (1:4)*5,fill=rainbow(4),cex = 0.5)
#summary(lm(log10(R2.alea.predict(N,Pmin,Qmax)/mean)[-1]~log10(Pmin)[-1],data=R2_obs.tab.PQ))
```
```{r}
with(data = R2_obs.tab.PQ,
plot(Pmin,
R2.alea.predict(N,Pmin,Qmax)-mean,
col=rainbow(4)[N/5],
pch=(1+(Qmax==Pmin)),
cex=0.5,
log="x")
)
abline(h=1)
abline(v=2)
abline(v=c(5,10,15,20)*4,lty=2,col=rainbow(4))
abline(v=c(5,10,15,20),lty=3,col=rainbow(4))
#abline(lm(log10(R2.alea.predict(N,Pmin,Qmax)/mean)[-1]~log10(Pmin)[-1],data=R2_obs.tab.PQ))
legend("topright",legend = (1:4)*5,fill=rainbow(4),cex = 0.5)
#summary(lm(log10(R2.alea.predict(N,Pmin,Qmax)/mean)[-1]~log10(Pmin)[-1],data=R2_obs.tab.PQ))
```
```{r}
with(data = R2_obs.tab.PQ,
plot(Pmin,
mean,
col=rainbow(4)[N/5],
pch=(1+(Qmax==Pmin)),
cex=0.5,
log="x")
)
abline(h=1)
legend("topright",legend = (1:4)*5,fill=rainbow(4),cex = 0.5)
```
=======
>>>>>>> ccb4d408c4936c553b8d097d94022285a7bc89a7
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment