異なる周期の三角関数の合成いろいろ

フーリエ変換の予備知識として
色々な三角関数の合成を可視化しておきます
あくまで準備段階ですので…

弧度法を度数法に変換します
x<-pi/180
y<-0:360
x<-x*y　

sin(x)、sin(2x) の合成
plot(sin(x)+sin(2*x),type="l",ylim = c(-2,2),ann = F)
lines(sin(x),type="l",col=2)　 #赤
lines(sin(2*x),type="l",col=3) #緑

sin(x)、sin(2x)、sin(3x)の合成
plot(sin(x)+sin(2*x)+sin(3*x),type="l",ylim = c(-2.5,2.5),ann = F)
lines(sin(x),type="l",col=2) 　#赤
lines(sin(2*x),type="l",col=3) 　#緑
lines(sin(3*x),type="l",col=4) 　#青

cos(x)、sin(3x)、-0.3cos(2x)の合成
plot(cos(x)+sin(3*x)-0.3*cos(2*x),type="l",ylim = c(-2.5,2.5),ann = F)
lines(cos(x),type="l",col=2) 　#赤
lines(sin(3*x),type="l",col=3) 　#緑
lines(-0.3*cos(2*x),type="l",col=4) 　#青

cos(nx)とsin(nx)の合成

グラフでイメージするa*cos(x)とb*sin(x)の合成
ベクトル(振幅)は√（a^2+b^2）
sun(nx)、cos(nx)を合成すると位相が変化するが、周期は変化なし

振幅　√(2^2+2^2)=2.828427 #sqrt(8)
1周期
plot(2*cos(y*x)+2*sin(y*x),type="l",ylim = c(-3,3),ann = F)
lines(2*cos(y*x),type="l",col=2)　#赤
lines(2*sin(y*x),type="l",col=3)　#緑

2周期
plot(2*cos(2*y*x)+2*sin(2*y*x),type="l",ylim = c(-3,3),ann = F)
lines(2*cos(2*y*x),type="l",col=2)　#赤
lines(2*sin(2*y*x),type="l",col=3)　#緑

　　　　　　　　　振幅は変化なし、周期のみが変化

振幅√(2^2+5^2)=5.385165 #sqrt(29)
2周期
plot(2*cos(2*y*x)+5*sin(2*y*x),type="l",ylim = c(-5.5,5.5),ann = F)
lines(2*cos(2*y*x),type="l",col=2)
lines(5*sin(2*y*x),type="l",col=3)

　　　　　　　　　　周期は変化なし、振幅のみが変化

三角関数のグラフ

角度の概念を復習します

弧度法
ラジアン(弧度)とは、「弧の長さ＝半径」となる角度を「1」とする角度の単位のこと。
単位記号は「rad」
R言語ではこの「弧度法」を使用しています．

度数法
円周を360等分したものを1とする角度の単位のこと。
単位記号は「°」

つまり2πr÷r が360度に値します
2πが360度

グラフを理解するために弧度法を度数法に変換します
x<-pi/180
90度はx*90 となります。

sin(90度)=sin(90*x)=sin(pi/2)=1

したがって0°から360°までは次のようになります
y=0:360　
x*y＝0°から360°まで

0°から360°までのsin波(黒),cos波(赤),sin波×cos波(緑),sin波×cos波(青)のグラフを重ねて書いてみます
plot(sin(y*x),type="l",ylim = c(-1.5,1.5),ann = F)　 　　#タイトルを全て消去 ann = F
lines(cos(y*x),type="l",col=2)
lines(cos(y*x)+sin(y*x),type="l",col=3)
lines(cos(y*x)*sin(y*x),type="l",col=4)
na<-c("sin波","cos波","sin波+cos波","sin波×cos波") 　#波に名前をつけます
legend(230,1.5,bty="n",na,text.col=1:4) 　　　　　　#bty="n" 凡例の枠無し
abline(h = 0,lty=2) 　　　　　　　　　　　　　　　 #ｙ軸のライン　(破線)

RでBland Altman Plotを描く

データ
国際医療福祉大学　下井研究室
信頼性の検討方法：相対信頼性と絶対信頼性
例1（http://shimoi.iuhw.ac.jp/reliability_3_systematicbias.html）
を使用させていただきました

x1
25.1 23.3 33.2 20.9 22.8 34.9 26.5 26.1 29.4 29.8 29.2 26.5 14.3 14.6 17.3 20.2 17.4 19.6 30.2 30.3 23.1 17.7 18.9 28.0 15.2 16.8 11.7 10.5 25.2 15.9

y1
17.6 21.5 25.6 10.5 16.4 31.2 20.8 27.0 32.8 30.9 22.9 24.0 18.2 14.9 13.2 17.2 19.5 15.8 21.2 26.1 14.1 14.8 12.3 17.6 14.3 8.2 11.3 11.8 23.4 10.6

まずはインストールを忘れずに

install.packages("BlandAltmanLeh")
library(BlandAltmanLeh)

ba <- bland.altman.stats(x1, y1)  
str(ba)　　　　　　　　　　　 #データを概観して必要な部分だけ抜き取ります
plot(ba$means, ba$diffs,ylim=c(-8,15))
abline(h=ba$mean.diffs, lty=2, col=2)
abline(h=ba$CI.lines[2], lty=2)　#誤差の許容範囲として最もあまい(optimisticな)範囲
abline(h=ba$CI.lines[5], lty=2)　#誤差の許容範囲として最もあまい(optimisticな)範囲

bland.altman.plot(x1,y1)

bland.altman.plot(x1,y1,conf.int = .95)

1標本ｔ検定のスクリプト

只今練習中～

2標本xとyについて
sumxy <- function(x,y,sxy)　　
{
xs<-summary(x)
ys<-summary(y)
xt<-t.test(x)
yt<-t.test(y)
print(xs);print(ys);print(xt);print(yt)
}
sumxy (x,y,sxy)
より

x<-c( );y<-c( )にデータを入れると、
各標本の要約・t値・p値・95%CIのみを出力させるプログラム

sumxy <- function(x,y,sxy)　　

{
cat1<-c("Min.","1st Qu","Median","Mean", "3rd Qu.","Max.")
xs<-summary(x)
ys<-summary(y)
cat(cat1,"\n","x",xs,"\n","y",ys,"\n\n",file = "TEST001") 　#各標本のまとめ

tx<-t.test(x)
cat2<-c("tvalue","pvalue","95%CI")
txt<- tx$statistic
txp<- tx$p.value
txc1<- tx$conf.int[1]
txc2<- tx$conf.int[2]
x1 <- sprintf("%.3f\t%.3f\t%.3f\t%.3f",txt,txp,txc1,txc2)

ty<-t.test(y)
tyt<- ty$statistic
typ<- ty$p.value
tyc1<- ty$conf.int[1]
tyc2<- ty$conf.int[2]
y1 <- sprintf("%.3f\t%.3f\t%.3f\t%.3f",tyt,typ,tyc1,tyc2)
cat(cat2,"\n","x",x1,"\n","y",y1,"\n",file = "TEST001", append =T) 　

}
sumxy (x,y,sxy)

繰り返し for

いつも忘れるので書いておこう

a = 0 に2を3回足す

a<-0

for(i in 1:3){
a<-a+2}

a

=6

#x=3に2を3回足す式
#3+2+2+2=9

myf <- function(x){
  x<-3　　 #ｘは3から開始
  for(i in 1:3){ 　　#x<-x+2でｘを書き直して３回繰り返し
    x<-x+2}
    return(x)
    }

myf()

#Xに2を3回加える式
#x+2+2+2
#xに5を代入すると5+2+2+2=11

myf<-function(x){
a<-x
for(i in 1:3){
a<-a+2}
return(a)
}
myf(5)

#Xに3を掛けることを3回繰り返す
#X*3*3*3
#xに5を代入すると
5*3*3*3

myf<-function(x){
for(i in 1:3){
x<-x*3}
return(x)
}
myf(5)

#10にｘを3回加える式
#10+x+x+x
#x=2のとき10+2+2+2=16

myf<-function(x){
a<-10
for(i in 1:3){
a<-a+x}
return(a)
}
myf(2)

相関係数の区間推定

サンプル：日本統計学会 (編集); 日本統計学会公式認定　統計検定２級対応 統計学基礎, 東京図書, 2012,p125

身長のデータ
大学生
x<-c(172,167,184,175,176,175,170,180,170,179,167,175,174,162,165,163,170,169,165,175)
父親
y<-c(165,165,178,176,150,171,172,175,170,156,163,170,165,160,163,170,163,165,160,172)

cor.test関数を使用します
同じ数の要素をもつベクトル２つの相関関係、また相関係数の 95% 信頼区間，p値を求めます．Pearsonの相関係数をもとめるのであれば次の式を実行するのみでOKです．

　　cor.test(x,y)

Pearson's product-moment correlation

data: x and y
t = 1.4587, df = 18, p-value = 0.1619
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1371017 0.6711054
sample estimates:
cor
0.3251458

cor.test(x,y)の引数
x, y：同じ大きさの数字ベクトル
alternative(hypothesis)： "two.sided", "greater", "less".
method ： "pearson(相関係数cor)", "kendall"(相関係数tau), "spearman"(相関係数rho).

　　　（r-de様のコメントを参考に修正しました）

次に上記関数結果の95％信頼区間[-0.1371017 0.6711054]の求め方を忘れないように残しておきます

大学生の標準偏差
Vx<-(sum( (x-mean(x) )^2) )/20
Sx<-sqrt(Vx)
Sx

父親の標準偏差
Vy<-(sum( (y-mean(y) )^2) )/20
Sy<-sqrt(Vy)
Sy

大学生と父親の共分散
Sxy<-(sum((x-mean(x) )*(y-mean(y) ) ) )/20

標本の相関係数r
Sxy/(Sx*Sy)
r=0.325

母相関係数ρの区間推定

標本の相関係数rの分布を考えてみます．
r＝0の場合
t=r*(sqrt(n-2)/sqrt(1-r^2) )は自由度n-2のt分布に従う

r≠0の場合（FisherのZ変換）
rの分布はρが0から離れると非対称性が強くなります．
Z={log((1+r)/(1-r) ) }/2と変換した場合、Zの分布は対称になります．
またZは近似的にN({log( (1+ρ)/(1-ρ) ) }/2,1/(n-3) )に従うことが知られています．

r=0.325をZ={log((1+r)/(1-r) ) }/2に代入してZ＝0.337を求めます．

{log((1+ρ)/(1-ρ) )＝Rとおくと

(0.337-R)＝±1.96×標準偏差
0.337－1.96／√(17)＜R＜0.337+1.96／√(17)
2×（0.337－1.96／√(17)）＜log((1+ρ)/(1-ρ)＜2×（0.337+1.96／√(17)）
‐0.138×2< log((1+ρ)/(1-ρ) )<0.671×2
したがってρの95％信頼区間は
(exp(-0.276)-1)/(exp(-0.276)+1)
-0.137
(exp(1.626)-1)/(exp(1.626)+1)
0.671

また勉強して更新していきますのでコメントお願いします．