プログラミング言語 R グラフ

画像のグラフを書きたいです。
R言語に詳しい方よろしくお願いします

No.5 ベストアンサー 回答者： kamiyasiro 回答日時： 2021/04/15 23:24

次に投稿された次元削減のご質問は、閉じられたのですね。

間に合わなくて残念です。
未完成ですが（図にも不一致がありますが）、スクリプトをこちらに載せておきます。ご参考まで。
他の回答者の方からの補足を頂きたいです。

install.packages()は現在コメントアウトしていますが、ライブラリが無い時は有効にしてインストールして下さい。

図は１個１個描いています。

ISOMAPはめちゃくちゃ時間がかかりますので、結果が出るのを根気よく待って下さい。

SIRは、目的変数をどうすれば良いかわからないので、テキトーです。

最後の２つのSIRは、方法を勉強しているところでした。ですから空欄です。

# ～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～

rm(list = ls())
par(ask = T)


# 2次元正規分布の分散共分散行列を与える

r <- 0 # 直交
Sigma <- matrix(c(1, r, r, 1), ncol = 2)

# 2次元正規分布の形状(高さ)を与える関数を定義する

#install.packages("mvtnorm")
library(mvtnorm)

Bivariate.Gauss <- function(x, y){
return(10 * dmvnorm(cbind(x, y), c(0, 0), Sigma, log=FALSE))
}

# 3次元の座標群を生成する

x <- y <- seq(-1.96, 1.96, length = 41) # 範囲と分割数でグラフは変わり得る
z <- outer(x, y, Bivariate.Gauss)
nrz <- nrow(z)
ncz <- ncol(z)

# パース図を描いて形状を可視化しチェックする～～～～～～～～～～～～～～～～～～

# Generate the desired number of colors from this palette
# レインボーカラーのうち、青⇔赤を使って色分けする

color <- rainbow(100)[70:1]

# Compute the z-value at the facet centres
# パース図の各ファセットの四隅の平均を計算してファセットのz軸値とする
# 具体的には、右下の点、左下の点、右上の点、左上の点を抽出して足す

zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]

# Recode facet z-values into color indices
# ファセットのz軸の値を離散化してレインボーカラーに対応させる

facetcol <- cut(zfacet, 70)

# パース図を描く

persp( x, y, z,
　theta = 30, phi = 30, # 回転させる角度
　expand = 1, # z軸の表示倍率(デフォルトは 1)
　col = color[facetcol], # 色の設定
　border = NA
)

# ～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～

# データをマージしておく

xyz <- data.frame(expand.grid(x, y), as.vector(z))
colnames(xyz) <- c("x", "y", "z")
xx <- as.matrix(scale(xyz))
color <- color[cut(as.vector(z), 70)]

# 次元縮約を掛ける前のx-y平面プロット

plot(xx[, 1:2], pch = 16, col = color,
xlim = c(-2, 2), ylim = c(-2, 2),
xlab = "", ylab = "")

# ～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～

# MDS

library(stats) # 組み込みライブラリ?

eurodist <- dist(xx)
result1 <- cmdscale(eurodist)
plot(result1[, 1:2], pch = 16, col = color,
xlim = c(-2, 2), ylim = c(-2, 2),
xlab = "", ylab = "", main = "MDS")


# PCA

result2 <- prcomp(xx)
plot(result2$x[, 1:2], pch = 16, col = color,
xlim = c(-2, 2), ylim = c(-2, 2),
xlab = "", ylab = "", main = "PCA")


# SIR

#install.packages("dr")
library(dr)

d <- mahalanobis(xyz, apply(xyz, 2, mean), cov(xyz))
data <- cbind(xyz, d)
result3 <- dr(d ~ ., data = data)
s <- update(result3, slice.function = dr.slices.arc)
plot(as.matrix(xyz) %*% s$evectors[,1:2], pch = 16, col = color,
xlim = c(-2, 2), ylim = c(-2, 2),
xlab = "", ylab = "", main = "SIR")


# ICA

#install.packages("fastICA")
library(fastICA)

result4 <- fastICA(xx, 2, alg.typ = "parallel", fun = "logcosh", alpha = 1,
method = "C", row.norm = FALSE, maxit = 200,
tol = 0.0001, verbose = TRUE)
plot(result4$S, pch = 16, col = color,
xlim = c(-2, 2), ylim = c(-2, 2),
xlab = "", ylab = "", main = "ICA")


# ISOMAP

#install.packages("vegan")
library(vegan)

eurodist <- dist(xx)
result5 <- isomap(eurodist, ndim = 2, k = 5)
plot(result5$points[, 1:2], pch = 16, col = color,
xlim = c(-3, 3), ylim = c(-3, 3),
xlab = "", ylab = "",, main = "ISOmap")


# K-PCA

#install.packages("kernlab")
library(kernlab)

result6 <- kpca(xx, kernel = "rbfdot", features = 2, kpar = list(sigma = 0.05))
plot(pcv(result6), pch=16, col = color,
xlim = c(-0.2, 0.2), ylim = c(-0.2, 0.2),
xlab = "", ylab = "", main = "Kernel PCA")


# K-SIR


# ISOSIR

この回答へのお礼

わざわざ毎回ありがとうございます。とても助かっています。
まだまだR言語の質問があると思いますので、お時間がありましたらこれからもお願いします。

お礼日時：2021/04/29 01:34

No.4 回答者： kamiyasiro 回答日時： 2021/04/15 16:24

すみません。

先日ご質問された中に、Rのカラーマップを描くって問題がありましたが、参考として挙げたHPのスクリプトの可読性が低かったので、分かりやすく書き直しました。（極力、変数を減らしました）
プロットの丸はキャラクタ番号21で描き、背景色を変更しています。
一度、動かしてみて下さい。

# R color spectrums

rm(list = ls())
n <- 24

main <- c("R color spectrums")
col.name <- c("rainbow", "heat.colors", "terrain.colors", "topo.colors", "cm.colors", "gray")
bg.col <- cbind(rainbow(n), heat.colors(n), terrain.colors(n), topo.colors(n), cm.colors(n), gray(1:n / n))

plot(0, 0, type = "n", axes = FALSE, xlab = "", ylab = "", main = main, xlim = c(0, n + 1), ylim = c(0, 6.5))

for(k in 1:6){
points(c(1:n), rep(k - 0.5, n), pch = 21, cex = 2.5, col = 8, bg = bg.col[, k])
text(n / 2, k, c(col.name[k]))
}

box()

No.3 回答者： kamiyasiro 回答日時： 2021/04/09 01:01

#2です。

この図の狙いは、「二山の離れ方」と「縦の破線の位置」を変更すると、統計的検定力（Statistical Power）＝１－β（水色部分）がどう変化するかということを示したいのだと思います。

そのためには、今のスクリプトでは、
・対立仮説Haの平均を帰無仮説H0の2.5σの位置にしていますが、それを変数に、
・有意水準（縦の破線）をH0の2σの位置（本来なら95％信頼限界は±1.96σ）にしていますが、それを変数にして、
それらを各々変化させたとき、αと(１－β)の関係がどうなるかということを調べるようにスクリプトを直さないといけません。
そういう課題なのではありませんか？

あと、図のタイトルの「Textbook-Style」が気になっているのですが・・・。そうじゃないスタイルがあるとでも言いたいのかしらん。
Texbookは両側検定（片側の有意水準＝α／2）だということでしょうかね。

No.2 回答者： kamiyasiro 回答日時： 2021/04/07 13:31

#1です。

注釈文字なしであれば、簡単に描けます。
高さとかの調整はご自分でお願いします。
分からなければ、「お礼」に書いてください。「補足」に書いても私のような投稿者には連絡が来ません。

rm(list = ls())

x <- seq(-4, 6.5, by = 0.05)

plot(0, 0, xlim = c(-4, 6.5), ylim = c(0, 0.5), axes = F,
xlab = "", ylab = "", pch = "")

#abline(h = 0, lty = 1)

polygon(c(seq(-4, 2, by = 0.05), 2), c(dnorm(seq(-4, 2, by = 0.05), 2.5, 1), 0), col = "#80400050", lty = 0)
polygon(c(2, seq(2, 6.5, by = 0.05)), c(0, dnorm(seq(2, 6.5, by = 0.05), 2.5, 1)), col = "skyblue", lty = 0)
polygon(c(2, seq(2, 6.5, by = 0.05)), c(0, dnorm(seq(2, 6.5, by = 0.05), 0, 1)), col = "steelblue", lty = 0)

lines(x, dnorm(x, 0, 1), col = 1)
lines(x, dnorm(x, 2.5, 1), col = "#804000")

abline(v = 2, lty = 2)

mtext("Statistical Power Plots, Textbook-style", line = 0)

No.1 回答者： kamiyasiro 回答日時： 2021/04/07 11:00

企業で統計を推進する立場の者です。

引き出し線を描いて注釈を入れるのは、手作業だと思います。
にも拘わらず、この図を描けという課題が出たのですか？