X
トップページ⇒数学
1002コメント395KB

【R言語】統計解析フリーソフトＲ 第6章【GNU R】 [無断転載禁止]©2ch.net

■ このスレッドは過去ログ倉庫に格納されています
0001１３２人目の素数さん
垢版 |
2017/08/03(木) 19:23:12.67ID:Hq1blL0O
R は統計計算とグラフィックスのための言語・環境です。
統計計算で重宝するデータ型や、複数要素を処理する演算や関数、
解析結果を表示するグラフィックなど、多彩な機能を提供します。

●関連サイト
The R Project
http://www.r-project.org/
RjpWiki
http://www.okada.jp.org/RWiki/
リンク集
http://www.okada.jp.org/RWiki/?%A5%EA%A5%F3%A5%AF%BD%B8
※前スレ
【R言語】統計解析フリーソフトＲ 第５章【GNU R】
http://rio2016.2ch.net/test/read.cgi/math/1380168442/
0780１３２人目の素数さん
垢版 |
2020/02/08(土) 19:51:16.44ID:+JttXwIS
実行したら全員が下船できるまで約1ヶ月になった。
コードにまだバグがあるかも。
0781１３２人目の素数さん
垢版 |
2020/02/09(日) 21:31:53.47ID:IDQ3+Iu7
Rのcsvの文字化け治らねえ
csvそのものをutf-8にしたら今度は不正なマルチバイトがあるってエラー
どうすればいいか教えて
0783１３２人目の素数さん
垢版 |
2020/02/09(日) 23:10:12.62ID:KodPd2ez
>>781
エスパーじゃないから、まず確認だけど

・Rを使ってる環境(OS)とRのバージョンは？
・CSVのコード変換に使ってるソフトは？
・CSVを読み込んでる関数と指定している引数は？
0785１３２人目の素数さん
垢版 |
2020/03/07(土) 09:46:33.43ID:3M7vmA2q
>>784
エラーがでた。R 3.6.3にアップデートしたけど
> library('remotes')
> remotes::install_github("GuangchuangYu/nCov2019", dependencies = TRUE)
Downloading GitHub repo GuangchuangYu/nCov2019@master
Skipping 1 packages not available: BiocStyle
Installing 5 packages: downloader, cowplot, maps, magick, BiocStyle
Error: Failed to install 'nCov2019' from GitHub:
(converted from warning) package ‘BiocStyle’ is not available (for R version 3.6.3)
0788１３２人目の素数さん
垢版 |
2020/03/07(土) 11:38:15.72ID:3M7vmA2q
We also developed a web app (http://www.bcloud.org/e/) with interactive plots and simple time-series forecasts.

という記述があるんだけど　 関数 nls による非線形回帰分析　で作成しているのかなぁ？
0789１３２人目の素数さん
垢版 |
2020/03/09(月) 20:05:18.91ID:0N1NTePA
面白い問題スレをRでシラミ潰しにカウントしてみた。

"縦n個、横n個のマス目のそれぞれに 1,2,3,...,n の数字を入れていく。
このマス目の横の並びを行といい、縦の並びを列という。
どの行にも、どの列にも、2つの対角線上にも同じ数字が1回しか現れない入れ方は何通りあるか求めよ。(2020京大文系　改)
n=4がオリジナル
"
library(gtools)
n=4
(a=permutations(n,n)) # ｎの順列
r=nrow(a)
f<-function(x){ # x=c(1,2,3,4) -> rbind(a[1],a[2],a[3],a[4])
ans=NULL
for(i in x){
ans=rbind(ans,a[i,])
}
return(ans)
}

b=permutations(r,n)
head(b) ; tail(b)
B=apply(b,1,f)
diag2 <- function(x){ # 行列xの右上から左下への対角線の配列を返す
n=nrow(x)
ans=numeric(n)
for(i in 1:n) ans[i]=x[i,n+1-i]
return(ans)
}
g <- function(v){ # matri(v,n)で列、対角線の要素がすべて異なればTRUEを返す
n=sqrt(length(v))
(x=matrix(v,n))
if(length(unique(diag (x))) < n) return(FALSE)
if(length(unique(diag2(x))) < n) return(FALSE)
flag=TRUE
for(i in 1:n){
if(length(unique(x[,i])) < n){
flag=FALSE
break
}
}
return(flag)
}
counter=NULL
for(i in 1:ncol(B)){
if(g(B[,i])) counter=append(counter,i)
}
length(counter)
counter
matrix(B[,counter[1]],n)
matrix(B[,counter[2]],n)
matrix(B[,counter[48]],n)
0790１３２人目の素数さん
垢版 |
2020/03/10(火) 07:11:26.36ID:H1fx2jVB
>>789
n=5 にするとメモリー不足でエラーがでた。
サンプリングでやってみる
# simulation
n=5
sim <- function(n=5){
v=as.vector(t(replicate(n,sample(n)))) # n × n行列要素をベクトル化
if(g(v,dia=TRUE)){ # diagonal latin squareなら
print(matrix(v,n)) # 行列化して表示
return(1)
}
return(0)
}
k=1e6
sum(replicate(k,sim(5)))
0791１３２人目の素数さん
垢版 |
2020/03/10(火) 12:07:20.32ID:KHQs8ybj
１億回モンテカルロしてようやく、1個みつかった。


[,1] [,2] [,3] [,4] [,5]
[1,] 2 5 4 3 1
[2,] 4 3 1 2 5
[3,] 1 2 5 4 3
[4,] 5 4 3 1 2
[5,] 3 1 2 5 4
0792１３２人目の素数さん
垢版 |
2020/03/11(水) 11:06:12.00ID:7KW6Eui5
嫌ね最近は総当たりアルゴリズムで枝刈り使えば実用的な時間で解決できるこ知らない人が多くて
0793１３２人目の素数さん
垢版 |
2020/03/15(日) 08:45:11.95ID:OTl1KJku
versionが3.6,xになって、str2lang とか　str2expressionという関数が追加されたんだな。

https://id.fnshr.info/2019/07/20/r3-6-0/
で知った

例

rep("1:6",5)
paste(rep("1:6",5),collapse=',')
str2lang(paste("expand.grid(",paste(rep("1:6",5),collapse=','),')'))
eval(str2lang(paste("expand.grid(",paste(rep("1:6",5),collapse=','),')')))

(s=rep("1:6",5))
(str=paste(s,collapse=','))
(lang=str2lang(paste("expand.grid(",str,')')))
eval(lang)
0794１３２人目の素数さん
垢版 |
2020/03/21(土) 17:23:30.68ID:RyI2Q/uv
The Clinical and Chest CT Features Associated with Severe and Critical COVID-19 Pneumonia

https://journals.lww.com/investigativeradiology/Abstract/publishahead/The_Clinical_and_Chest_CT_Features_Associated_with.98832.aspx

の付表にこんなデータがあった

# Parameter Total (n=83) Severe/critica l Group (n=25) Ordinary Group (n=58) P Value
# Number of lobes Median (interquartile range, IQR)
# 5 (4-5) 5 (5-5) 5 (3-5) 0.003

重症群２５例　中等症群５８例でCTで肺炎像が確認できた肺葉数(１〜５）が
全体で中央値が５　第一四分位数(少ない方から数えて25%の値）が４、第三四分位数(少ない方から数えて75%の値）が５
重症群では中央値５、第一四分位数５、第三四分位数５
中等症群では中央値５、第一四分位数３、第三四分位数５で
p値が０．００３であるという。
０．００２５〜０．００３５までとして、

この条件に見合う８３例をリストアップすると
> re[100,][1:25]
[1] 1 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
> re[100,][26:83]
[1] 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
[48] 5 5 5 5 5 5 5 5 5 5 5
という類の組み合わせで
19504通りと計算された。

重症群の肺葉数の平均は４．７６
中等症群の平均は４．０４
と計算された。

タイ値がこれだけ多いデータにWilcox検定を使う意味がわからんな。
0795１３２人目の素数さん
垢版 |
2020/03/22(日) 08:57:06.28ID:liILqu/N
# 感度0.7　特異度0.9のPCR検査で1000人中10人が陽性だったときの有病率の信頼区間をシミュレーションして計算するスクリプト

# アルゴリズムには自信がないｗ

sim <- function(n=1000,r=10,sen=0.7,spc=0.9,k=1e6,print=TRUE){
# n:検査件数 r:検査陽性数 sen=感度 spc=特異度 k:発生乱数の数,print:グラフ描画
prev=rbeta(k,1+r,1+n-r) # 有病率の事前分布を一様分布と仮定

PPV=prev*sen/( prev*sen + (1-prev)*(1-spc) ) # 検査特性を加味してPPV計算
prev.post = PPV*(1-spc)/(sen - PPV*(sen+spc-1)) # そのPPVから有病率を計算
# NPVでも同じ結果
# NPV=(1-prev)*spc/( (1-prev)*spc + prev*(1-sen)) # 検査特性を加味してNPV計算
　# prev.post = (1-NPV)*spc/(spc - NPV*(sen+spc-1)) # そのNPVから有病率を計算
mean=mean(prev.post) # 期待値
median=median(prev.post) # 中央値
mode=density(prev.post)$x[which.max(density(prev.post)$y)] # 最頻値
LU=unlist(HDInterval::hdi(prev.post))[1:2] # 95%信頼区間
re=c(mean=mean,median=median,mode=mode,LU) # 事後有病率の代表値
if(print){ # ヒストグラム描画
par(mfrow=c(2,2))
hist(prev,freq=F,xlim=c(0,1)) ; lines(density(prev),lwd=2)
hist(PPV,freq=F,xlim=c(0,1)) ; lines(density(PPV),lwd=2)
hist(prev.post,freq=F,xlim=c(0,1)) ; lines(density(prev.post),lwd=2)
BEST::plotPost(prev.post)
par(mfrow=c(1,1))
}
print(round(re,4)) # 概算値表示
invisible(list(re,prev.post)) # 代表値と事後有病率を非表示で返す
}
sim(100,1)
sim(1000,10)
sim(10000,100)
sim(100000,1000)
0797１３２人目の素数さん
垢版 |
2020/03/23(月) 14:55:08.54ID:iIfIhG5+
MCMCで解決してみた。

【富山県最強伝説】新型コロナウイルスPCR検査件数　54人　陽性0人

ある集団から54人を無作為に選んでPCR検査したら陽性0であった。
PCR検査の感度0.7　特異度0.9としてこの集団の有病率の期待値と95%信頼区間を推測せよ。


#----------------------------------------------------------------------------
# 感度SEN, 特異度SPCの検査でN人中X人が陽性であったときの推定有病率prevalence
# 弱情報事前分布はprevalence ~ dunif(0,UL) , UL:上限
#----------------------------------------------------------------------------

PCRj <- function(N,X,UL=1,SEN=0.7,SPC=0.9,verbose=TRUE){ # UL:upper limit of dunif(0,UL)
library(rjags)
modelstring=paste0('
model
{
x ~ dbin(p,n)
p <- prev*sen + (1-prev)*(1-spc)
prev ~ dunif(0,',UL,')
}
')
if(verbose & UL!=1) cat(modelstring)
writeLines(modelstring,'TEMPmodel.txt')
dataList=list(n=N,x=X,sen=SEN,spc=SPC)
jagsModel = jags.model( file="TEMPmodel.txt" ,data=dataList,quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples( jagsModel ,
variable=c("prev","p"), n.iter=1e6, thin=10)
js=as.matrix(codaSamples)
BEST::plotPost(js[,'prev'],xlab='prevalence',showMode = TRUE) ; lines(density(js[,'prev']),col='skyblue')
round(c(mean=mean(js[,'prev']),HDInterval::hdi(js[,'prev'])[1:2]),10)
}

PCRj(100,1)
PCRj(1000,10)
PCRj(10000,100)
PCRj(100000,1000)

PCRj(54,0)
PCRj(54,0,UL=0.1)
0798１３２人目の素数さん
垢版 |
2020/04/03(金) 19:48:06.63ID:MNrUuJMd
a b c d e f g h i j
1 -0.0377710 0.02994835 0.1779426 0.29057472 0.13252796 0.3279709 -0.3534731 -0.54812874 -0.3869768 -0.1233967
2 0.2707033 0.09238235 -0.1042958 -0.08485129 -0.02590955 0.2159095 -0.3338677 0.00865952 -0.1575856 0.0953795
0799１３２人目の素数さん
垢版 |
2020/04/03(金) 19:50:07.75ID:MNrUuJMd
上のDFを

library(psych)
ICC(DF）

で級内相間求めようとすると変な値でるんだけどなんででしょう？
0800１３２人目の素数さん
垢版 |
2020/04/03(金) 22:01:48.08ID:1EdxzcOL
> DF
V1 V2
a -0.03777100 0.27070330
b 0.02994835 0.09238235
c 0.17794260 -0.10429580
d 0.29057472 -0.08485129
e 0.13252796 -0.02590955
f 0.32797090 0.21590950
g -0.35347310 -0.33386770
h -0.54812874 0.00865952
i -0.38697680 -0.15758560
j -0.12339670 0.09537950
> library(psych)
> ICC(x=DF,lmer=FALSE)
Call: ICC(x = DF, lmer = FALSE)

Intraclass correlation coefficients
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.36 2.1 9 10 0.13 -0.18 0.74
Single_random_raters ICC2 0.34 2.0 9 9 0.17 -0.26 0.74
Single_fixed_raters ICC3 0.32 2.0 9 9 0.17 -0.24 0.72
Average_raters_absolute ICC1k 0.53 2.1 9 10 0.13 -0.43 0.85
Average_random_raters ICC2k 0.51 2.0 9 9 0.17 -0.69 0.85
Average_fixed_raters ICC3k 0.49 2.0 9 9 0.17 -0.63 0.84

Number of subjects = 10 Number of Judges = 2
0801１３２人目の素数さん
垢版 |
2020/04/03(金) 23:40:27.04ID:MNrUuJMd
>>800
すいません
書き方が良くなかったです。
subject が１，２
judge a,b,c,dをきちんとfirst,second,third,fourth
と書くべきでした。
subjectとjudgeはあってると思うのです。
0802１３２人目の素数さん
垢版 |
2020/04/03(金) 23:42:39.27ID:MNrUuJMd
例えば
first second third fourth
0 0.05 0.1 0.15
0 0 0.05 0.05
0.15 0.1 0.05 0.1
0 0.05 0.2 0.05
0.1 0.15 0.1 0.05
0.15 0.05 0.1 0.05

というデータで

Call: ICC(x = over_sight)

Intraclass correlation coefficients
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.0068 1 5 18 0.43 -0.19 0.48
Single_random_raters ICC2 0.0068 1 5 15 0.44 -0.19 0.48
Single_fixed_raters ICC3 0.0068 1 5 15 0.44 -0.19 0.48
Average_raters_absolute ICC1k 0.0268 1 5 18 0.43 -1.70 0.79
Average_random_raters ICC2k 0.0268 1 5 15 0.44 -1.82 0.79
Average_fixed_raters ICC3k 0.0268 1 5 15 0.44 -1.82 0.79

Number of subjects = 6 Number of Judges = 4
0803１３２人目の素数さん
垢版 |
2020/04/03(金) 23:44:30.81ID:MNrUuJMd
となるのですが、ICC1が0.0068ってこんなに低くなるのでしょうか？
結果のとり得る値が1〜0なので最大15％ぐらいの変動幅で、
制度としてはかなりいい検査方法だとおもったのですが・・・
0804１３２人目の素数さん
垢版 |
2020/04/04(土) 00:07:42.72ID:npCU99bV
少し理解できました。
範囲制約性の問題があるのですね。
今回のサンプルがみんな低い得点で分散が低いからICCが低く算出されていました。
ダミーデータ
1,1,1,1
を追加したらICC=0.99となりました。
得点の高いサンプルを追加しないとICCで信頼性は測れないのですね。
0805１３２人目の素数さん
垢版 |
2020/04/04(土) 00:09:17.81ID:npCU99bV
いま有病者と健常者の両方を計測する検査法で健常者のみのデータでICCを算出しようとしているのが間違いということですね。

上記のような場合で、健常者のみでも信頼性を算出する方法っていうのはないのでしょうか？
0806１３２人目の素数さん
垢版 |
2020/04/04(土) 21:44:33.98ID:ZFu90Xbq
>>804
DF=rbind(
c(0,0.05,0.1,0.15),
c(0,0,0.05,0.05),
c(0.15,0.1,0.05,0.1),
c(0,0.05,0.2,0.05),
c(0.1,0.15,0.1,0.05),
c(0.15,0.05,0.1,0.05),
c(0.0001,0.0001,0.0001,0.0001))

とダミーに0.0001を加えても
> DF2ICC1(DF)
ICC(1,1) ICC(1,k)
0.2462707 0.5665262

> psych::ICC(DF)
type ICC F df1 df2 p lower bound
Single_raters_absolute ICC1 0.25 2.3 6 21 0.072 -0.027
Single_random_raters ICC2 0.25 2.3 6 18 0.079 -0.028
Single_fixed_raters ICC3 0.25 2.3 6 18 0.079 -0.034
Average_raters_absolute ICC1k 0.57 2.3 6 21 0.072 -0.115
Average_random_raters ICC2k 0.57 2.3 6 18 0.079 -0.122
Average_fixed_raters ICC3k 0.57 2.3 6 18 0.079 -0.154

となるから、高値データでなくてもいいみたい。
0810１３２人目の素数さん
垢版 |
2020/04/25(土) 07:29:33.43ID:R/UD6QQG
>>809
これやってからでないとパッケージのコンパイルができない。

After installation is complete, you need to perform one more step to be able to compile R packages:
you need to put the location of the Rtools make utilities (bash, make, etc) on the PATH.
The easiest way to do so is create a text file .Renviron in your Documents folder which contains the following line:

PATH="${RTOOLS40_HOME}\usr\bin;${PATH}"

You can do this with a text editor, or you can even do it from R like so:

writeLines('PATH="${RTOOLS40_HOME}\\usr\\bin;${PATH}"', con = "~/.Renviron")
0811１３２人目の素数さん
垢版 |
2020/04/25(土) 08:56:31.54ID:R/UD6QQG
>>809
libraryのフォルダで　dir /B とやって既にinstall済のリストを作って
install.packages("rstan","ggplot2",....)とやった方がいいな。
5000ほどパッケージがあるらしいけど、全部インストールしたらどれくらいの容量が必要なんだろう？
0814１３２人目の素数さん
垢版 |
2020/04/26(日) 22:22:15.60ID:tZeixXMR
高血圧で不安だという老人の経過観察入院1件。発熱してなくてほっとした。
ステーキハウスがテイクアウトを始めたからインセンティブはその足しにしようっと。
0815１３２人目の素数さん
垢版 |
2020/05/03(日) 15:15:43.67ID:ICXX0aJG
>>812
コンパイルなしでやってみた。

サイズ：19.6 GB (21,082,530,395 バイト)

ディスク上のサイズ：20.3 GB (21,805,858,816 バイト)

ファイル数: 527,982、フォルダー数: 118,521
0816１３２人目の素数さん
垢版 |
2020/05/03(日) 15:37:44.00ID:f6MQMkc5
Ubuntu 20.04でRcmdrがビルドできないんですけど、これはソフトウェアが20.04に対応するまで待つしかないですか？
0819１３２人目の素数さん
垢版 |
2020/05/07(木) 01:08:38.98ID:VnQvkZ57
4.0になってから
col=番号指定での色が目に優しくなったみたい。

hist(runif(100),col=2)
hist(runif(100),col=rgb(1,0,0,1))
hist(runif(100),col='red')
0823１３２人目の素数さん
垢版 |
2020/05/07(木) 23:39:55.93ID:8IDmAqjy
リリースノートより

The palette() function has a new default set of colours (which are less saturated and have better accessibility properties).
There are also some new built-in palettes, which are listed by the new palette.pals() function.
These include the old default palette under the name "R3". Finally, the new palette.colors() function allows a subset of colours to be selected from any of the built-in palettes.
0825１３２人目の素数さん
垢版 |
2020/05/08(金) 00:23:14.39ID:CfFk/Uw1
binomにある信頼区間をグラフ表示させてみた。

binom.ci <- function(x,n,cl=0.95){
PEci=binom::binom.confint(x,n,conf.level=cl)
y=nrow(PEci):1
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
par(mar=c(3,3,2,2))
par(oma=c(0,4,0,0))
plot(PEci$upper,y,xlim=c(min(PEci$lower),max(PEci$upper)),type="n",
yaxt="n",ylab="",xlab="Confidence Interval",
main=paste("C.I. for", x,"out of",n))
points(PEci$upper,y,pch="|") ; points(PEci$lower,y,pch="|")
segments(PEci$lower,y,PEci$upper,y,lwd=3,
col=sample(colours(),1))
abline(v=PEci$mean[1],lty=3)
axis(2, at=y, labels=PEci$method,as.vector,las=2,cex=0.75)
PEci
}
binom.ci(3,312)
0826１３２人目の素数さん
垢版 |
2020/05/08(金) 00:48:23.14ID:CfFk/Uw1
>>824
情報、ありがとうございます。

?paletteのヘルプにあるスクリプトを実行すると体感できました。

## Demonstrate the colors 1:8 in different palettes using a custom matplot()
sinplot <- function(main=NULL) {
x <- outer(
seq(-pi, pi, length.out = 50),
seq(0, pi, length.out = 8),
function(x, y) sin(x - y)
)
matplot(x, type = "l", lwd = 4, lty = 1, col = 1:8, ylab = "", main=main)
}
sinplot("default palette")

palette("R3"); sinplot("R3")
palette("Okabe-Ito"); sinplot("Okabe-Ito")
palette("Tableau") ; sinplot("Tableau")
palette("default") # reset
0827１３２人目の素数さん
垢版 |
2020/05/09(土) 20:51:55.34ID:oDT7bFgO
rjagsのコードに

theta[i,j] <- min(1, exp(-alpha[i]*t[j]) + beta[i])といいれて、確率を１以下にして二項分布させてようとすると

Node inconsistent with parents　というエラーがでた。　stackoverflow.comも明確な答が見いだせず、あれこれいじって

どうも確率1や0で二項分布する可能性があるとエラーがでるみたい、というのを発見。

1を0.99999とするとエラーが回避できた。

1を　1 - 1e-16としてもエラー回避できた。
0828１３２人目の素数さん
垢版 |
2020/05/18(月) 08:52:51.27ID:s/hrJxRo
mathematicaは内部どうなってんだか知らないが第何種ベッセル関数とか入ってくるめちゃくちゃ重い積分を計算してくれた
0829１３２人目の素数さん
垢版 |
2020/05/25(月) 12:44:17.51ID:3k8+9x6G
"
ある人物Dが新型コロナ肺炎に罹患したとする。
行動調査によって発症前にキャバクラに行っており
接客したキャバ嬢Cが人物D発症の２日後に発症していたことがわかった。
キャバ嬢Cは人物Dから移されたと主張して１億円の賠償を求めている。
潜伏期間には幅がありキャバ嬢Cから移された可能性もあると主張して
その確率を計算して賠償金を値切りたい。
いくら値切れるか計算せよ。
"
潜伏期間が従う分布はこれ
#--- incubation period ---
# from Li et al NEJM 2020
# lognormal mean = 5.2
ln_par1 = 1.434065
ln_par2 = 0.6612

Gt <- function(delay){
C=rlnorm(1e6,ln_par1,ln_par2)
D=rlnorm(1e6,ln_par1,ln_par2)
mean(C-D > delay)
}
Gt(2)
0830１３２人目の素数さん
垢版 |
2020/05/25(月) 12:48:45.94ID:3k8+9x6G
怖い本で紹介されていたパッケージ　polspline　は秀逸。
発生させた乱数から確率密度関数を作ってくれる。
複雑な分布だと収束できず作れないこともあるけど

対数正規分布の差の分布を乱数発生させてpdfが作れた。

library(polspline)
delay=2
ln_par1 = 1.434065
ln_par2 = 0.6612
c=rlnorm(1e5,ln_par1,ln_par2)
d=rlnorm(1e5,ln_par1,ln_par2)
hist(c-d, freq=F,col='white',breaks=100,ylim=c(0,0.11),
xlim=c(-30,30),ann=F,axes=F) ; axis(1)
fit=logspline(c-d)
curve(dlogspline(x,fit),-30,30,ann=F,bty='n',add=T)
segments(delay,0,delay,dlogspline(delay,fit),pch=19,col=4)
curve(dlogspline(x,fit),delay,30,add=T,type='h',col=4)
curve(plogspline(x,fit),-30,30,lty=3,bty='n')
1-plogspline(delay,fit)
fn <- function(delay) 1- plogspline(delay,fit)
curve(fn(x),0,14, bty='n' ,xlab='Delay', ylab='Probability')
0832１３２人目の素数さん
垢版 |
2020/07/25(土) 08:42:04.47ID:1FL39oBq
Rのグラフの出力をGIF動画にしたかったので検索したら、このページが役立った。

http://cse.naro.affrc.go.jp/takezawa/r-tips/r/36.html

mytest <- function() {
ANS <- readline("1+2? ")
if (ANS == "3") cat("Correct!\n")
else cat("Wrong...\n")
}
mytest()


mymessage <- function() {
Sys.sleep(3)
cat("Hello.\n")
}
mymessage()

尚、GIF作製には　
https://www.screentogif.com/
がお手軽だった。
0833１３２人目の素数さん
垢版 |
2020/07/26(日) 10:34:51.02ID:FqjSGBlb
他人のコードを読んでいたら、

'%&%' <- function(x,y) paste0(x,y)

という関数が定義されていた。

これは便利と思った。

> '√2 =' %&% sqrt(2) %&% ', π = ' %&% pi
[1] "√2 =1.4142135623731, π = 3.14159265358979"
0837１３２人目の素数さん
垢版 |
2020/07/29(水) 12:22:51.45ID:3HY+hEHZ
面白いけど使いどころが難しそう
%>%パイプとの親和性が無いのが痛い
0840１３２人目の素数さん
垢版 |
2020/07/30(木) 18:02:20.97ID:8iPZt0Nw
>>833
それって意味があるの？
> paste0('√2 =', sqrt(2), ', π = ', pi)
[1] "√2 =1.4142135623731, π = 3.14159265358979"
普通にpaste0のまま使った方がタイプ量が圧倒的に少ないが。
0841１３２人目の素数さん
垢版 |
2020/07/30(木) 20:10:10.19ID:417El1m4
>>840
括弧の対応で混乱しそうなときにちょっと役にたつかも。

こんな感じ

eval(str2lang("expand.grid(" %&% paste0(rep("1:2",6),collapse=',') %&% ")"))
0842１３２人目の素数さん
垢版 |
2020/07/31(金) 07:51:08.50ID:0a0zx1wG
>>841
うーん、その例だとイマイチに思える

rerun(6, 1:2) %>% expand.grid
で済むものをわざわざ複雑にしているだけに感じる
0843１３２人目の素数さん
垢版 |
2020/07/31(金) 13:51:36.15ID:ImynyFWQ
>>842
無理やり、例示した感は否めないな。
標準装備でなら、
expand.grid(replicate(6,1:2,simplify = FALSE))
ですむし。
0844１３２人目の素数さん
垢版 |
2020/07/31(金) 23:36:09.03ID:0a0zx1wG
paste()やstr_c()以外を使うなら、sprintf()かglue::glue()が候補か
特にカンマが文字列に含まれるときにはpaste()より見やすくなるはず
sprintf("√2 = %s, π = %s", sqrt(2), pi)
glue("√2 = {sqrt(2)}, π = {pi}")
0845１３２人目の素数さん
垢版 |
2020/08/01(土) 18:51:48.81ID:2TFdzAyg
>>844
確かに便利ですね。
カンマを　','で入力するのは混乱のもとでしたから。
有用な情報、ありがとうございます。'

> glue::glue("√2 = {sqrt(2)}, sin(π/3) = {sin(pi/3)}")
√2 = 1.4142135623731, sin(π/3) = 0.866025403784439
0846１３２人目の素数さん
垢版 |
2020/08/06(木) 00:01:24.75ID:rivMs4GJ
初心者質問で申し訳ありません。
ある関数を作りたいのですが、その関数に入れる数値の数はランダムな場合はどのように書けばいいのでしょうか。
例えば、入れた要素をすべて掛け合わせる関数が作りたい場合、そこに入れる要素が（２，３，４）でも、（２，３，４，５，６）でも反応するような関数の書き方を教えてくださると助かります。
0848１３２人目の素数さん
垢版 |
2020/08/06(木) 01:18:04.99ID:rivMs4GJ
874さん。
847さん、回答ありがとうございます。
確かにベクトルで何かするのはわかるのですが、その部分がわからないのです。
出来れば下に問題を書くので、コードを書いていただけると非常に助かります。

問、正の整数であるベクトルを入力すると、そのベクトルの要素をすべて掛け合わせた値を結果として返す関数を作りなさい。
0849１３２人目の素数さん
垢版 |
2020/08/06(木) 02:03:47.23ID:V+bu3PeG
>>848
問題文これだけ？他に条件はないの？
標準の関数にあるけどいいのかな？

関数fを作るとして、一番いいのは標準の関数に余計な変更を加えないことだから、これが答えだけど…
f <- prod
0850１３２人目の素数さん
垢版 |
2020/08/06(木) 02:11:09.57ID:rivMs4GJ
849さん、返信ありがとうございます。

実は問２がありまして、たぶんそれは既存の関数にはないと思うので、そちらも教えていただけると非常に助かります。

問２、要素がすべて整数であるベクトルを入力すると、そのベクトルの要素の絶対値をすべて掛け合わせた値を結果として返す関数を作成せよ。
0854１３２人目の素数さん
垢版 |
2020/08/06(木) 05:32:55.68ID:meNNWVIo
>>850
f <- function(...){
v=c(...)
ans=1
for(i in v) ans=ans*i
ans=ifelse(ans>0,ans,-ans)
ans
}

> f(-1,2,3)
[1] 6
> f(1 -2,-3,-4,-5)
[1] 60
0855１３２人目の素数さん
垢版 |
2020/08/06(木) 05:34:50.08ID:meNNWVIo
>>850
入力はベクトルだったから、

> f <- function(v){
+ ans=1
+ for(i in v) ans=ans*i
+ ans=ifelse(ans>0,ans,-ans)
+ ans
+ }
> f(c(-1,2,3))
[1] 6
> f(c(1 -2,-3,-4,-5))
[1] 60
0856１３２人目の素数さん
垢版 |
2020/08/06(木) 06:04:14.18ID:meNNWVIo
再帰だとこんな感じかな？

f <- function(...){
v=c(...)
n=length(v)
sub<-function(n){
if(n==0) return(1)
else v[n]*Recall(n-1)
}
sub(n)
}
> f(2,3,4)
[1] 24
> f(2,3,4,5,6)
[1] 720
0859１３２人目の素数さん
垢版 |
2020/08/06(木) 08:49:09.44ID:E/ufUouA
再帰版はこんな感じかな。
rec.prod <- function(x) {
# 再帰により x の要素の絶対値の積を求める。
# 実際には prod(abs(x)) が良い。
l <- length(x) # ベクトルの要素数
if (l == 0) # 空のベクトルなら
return(1) # 積は1
x1.a <- abs(x[1])
if (l == 1)
x1.a
else
x1.a * rec.prod(x[-1])
}
0860１３２人目の素数さん
垢版 |
2020/08/06(木) 11:42:21.63ID:rivMs4GJ
質問者本人です。
皆様のお力添えのおかげで無事問題を解くことができました。
初心者質問にやさしく回答して頂いて感謝しかありません。

本当にありがとうございました。
0861１３２人目の素数さん
垢版 |
2020/08/06(木) 19:28:54.17ID:meNNWVIo
>>859
この方が>856みたいに関数の中で関数を定義とかなくてカッコいいと思ったのでsystem.timeで時間を比較しようと思って実験を始めたら、
何故か、rec.prodの方がoverflowしやすいなぁ。
頭からかけても尻からかけても同じじゃないかと思うんだが。

> prod(1:100)
[1] 9.332622e+157

> rec.prod <- function(x) {
+ # 再帰により x の要素の絶対値の積を求める。
+ # 実際には prod(abs(x)) が良い。
+ l <- length(x) # ベクトルの要素数
+ if (l == 0) # 空のベクトルなら
+ return(1) # 積は1
+ x1.a <- abs(x[1])
+ if (l == 1)
+ x1.a
+ else
+ x1.a * rec.prod(x[-1])
+ }
> rec.prod(1:100)
[1] NA
Warning message:
In x1.a * rec.prod(x[-1]) : NAs produced by integer overflow

> f <- function(...){
+ v=c(...)
+ n=length(v)
+ sub<-function(n){
+ if(n==0) return(1)
+ else v[n]*Recall(n-1)
+ }
+ sub(n)
+ }
> f(1:100)
[1] 9.332622e+157
0862１３２人目の素数さん
垢版 |
2020/08/06(木) 20:06:06.68ID:meNNWVIo
比較のために、関数名の長さを同じにして、時間のかかるRecallを関数名にしてsystem.timeで比較してみた。

> rm(list=ls())
> rp <- function(x) {
+ # 再帰により x の要素の絶対値の積を求める。
+ # 実際には prod(abs(x)) が良い。
+ l <- length(x) # ベクトルの要素数
+ if (l == 0) # 空のベクトルなら
+ return(1) # 積は1
+ x1.a <- abs(x[1])
+ if (l == 1)
+ x1.a
+ else
+ x1.a * rp(x[-1])
+ }
> f1 <- function(...){
+ v=c(...)
+ n=length(v)
+ sub<-function(n){
+ if(n==0) return(1)
+ else v[n]*Recall(n-1)
+ }
+ sub(n)
+ }
> f2 <- function(...){
+ v=c(...)
+ n=length(v)
+ sub<-function(n){
+ if(n==0) return(1)
+ else v[n]*sub(n-1)
+ }
+ sub(n)
+ }


> system.time(replicate(1e5,prod(1:10)))
user system elapsed
0.11 0.00 0.11
> system.time(replicate(1e5,rp(1:10)))
user system elapsed
2.64 0.00 2.64
> system.time(replicate(1e5,f1(1:10)))
user system elapsed
2.44 0.01 2.47
> system.time(replicate(1e5,f2(1:10)))
user system elapsed
1.86 0.00 1.87
0863１３２人目の素数さん
垢版 |
2020/08/06(木) 20:31:17.99ID:E/ufUouA
それはね、856は厳密な意味での再帰ではないから速いしリ、ソースを食わないんだよ。ループを再帰風に書いただけだから。それよりもループのほうが速いに決まっている。
この問題を再帰で解くのはムダ以外の何ものでもない。
0864１３２人目の素数さん
垢版 |
2020/08/06(木) 20:33:25.51ID:w8yiRRv5
>>863
喩えれば、
親分が子分を再帰で酷使しているだけだから親分は消耗しない
という理解でいい？
0865１３２人目の素数さん
垢版 |
2020/08/06(木) 20:47:18.11ID:E/ufUouA
全然違う。
再帰は、複雑な問題をより小さい問題に分割することがキモだ。
856は元のベクトルはそのまま変わらずに存在していて、ある特定の要素を一つずつ処理しているからループとやってることは変わりない。
859は、一つ短いベクトルを自分自身に処理させているから真の再帰なのだ。
0866１３２人目の素数さん
垢版 |
2020/08/06(木) 22:25:29.65ID:meNNWVIo
>>865
f自身は再帰関数じゃないけど、呼び出しているsubは再帰関数では？

こういうのと同じ。
f <- function(x){　# 数列の和の階乗を返す
m=sum(x)
sub <- function(n){
if(n==1) return(1)
else{
return(n*Recall(n-1))
}
}
fact(m)
}

> f(c(1,2,3))
[1] 720
0867１３２人目の素数さん
垢版 |
2020/08/06(木) 22:50:49.72ID:XP8VEZbb
>>865
(ご入力訂正)
f自身は再帰関数じゃないけど、呼び出しているsubは再帰関数では？

こういうのと同じ。
f <- function(x){　# 数列の和の階乗を返す
m=sum(x)
sub <- function(n){
if(n==1) return(1)
else{
return(n*Recall(n-1))
}
}
sub(m)
}

> f(c(1,2,3))
[1] 720
0868１３２人目の素数さん
垢版 |
2020/08/07(金) 21:08:59.64ID:S7qsZ31k
856はsub()の引数がvでなくnなことと、nの入力に対して関数外のv[n]が帰ってくるのが、少し気になるかも
856が再帰っぽくないってのは、859はxを分割・再帰していってるのに対して、856はnを減じていってその都度vのn番目の値を参照してるところかな

それから、859がrec.prod(1:13)であっさりオーバーフローしてしまうのは、integer型だから
あと実行速度が少し遅いのは、xl.aへの代入操作と、if(l ==1){}の余分な条件分岐のせいだね

その辺を修正してやれば、こんな感じかな
# 再帰で問２
f <- function(x){
if(length(x) == 0) return(1)
abs(x[1]) * f(x[-1])
}
0869１３２人目の素数さん
垢版 |
2020/08/08(土) 05:40:19.06ID:2ggSSq05
>>868
勉強になりました。

ついでにwhile版を書いてみました。
f <- function(x){
ans=1
while(length(x)){
ans=ans*abs(x[1])
x=x[-1]
}
ans
}
0870１３２人目の素数さん
垢版 |
2020/08/08(土) 09:18:18.60ID:FvlAeRBY
>>868
859が遅い根本の原因は、x[-1] が生成されること。
他の指摘は、細かい高速化の工夫だから初心者に薦めるのはどうかと。
現に869で全く実用的に意味のないコードを書いているし。
題意を素直に表したほうがよいのでは？
絶対値の積は、積の絶対値と等しいから、absを呼ぶのは1回で済むが、何も説明無しにそうするのは私的にはアウトだ。メンテナンスの時に自分や他人が苦しむから。
0871１３２人目の素数さん
垢版 |
2020/08/08(土) 11:01:30.91ID:559rm1Tc
>>870
意味不明なこと言いだしてるが大丈夫か、落ち着け
859は、無駄な条件分岐あるし、すぐオーバーフローするし、コードが冗長だしで、それこそ初心者以外にもNGだぞ

絶対値の積うんぬんに関しては、859と同じでabsを毎回呼んでる、余計な代入操作が無いだけだ
落ち着いてコードを見ろ

速度うんぬんの話はこれまで話題に出てきたように、856と比べての話だ、話についてこれないなら無理に入ってこなくていい
別に高速化を目的としてるわけじゃなく、859から冗長な部分やオーバーフローの原因を取り除いた結果高速化されたってだけだ

だが、速度が遅い原因がxl.aへの代入でなく、x[-1]作成が原因てのはそのとおりだ
文章を書いてる途中で、オーバーフローの原因（代入操作はこっちだ）と勘違いしてしまったようだ
これは訂正に感謝だ
0872１３２人目の素数さん
垢版 |
2020/08/08(土) 11:24:32.96ID:FvlAeRBY
>>871
おまえこそもちつけ。
再帰は、再帰の停止条件と、自分自身が行う作業に分離できる。
この場合の停止条件は、要素が1つになったときと考えるのが自然だろう。
数学的には要素が0でもよいし、Rの場合はそれでうまくいくから問題はない。だが、そういった詳細を何も説明せずにこれでよい、と示すのはどうかな？
あと、絶対値を取ったのを代入しているのは、この作業が複雑になった場合の修正の確実さを優先するためで、優秀なバイトコンパイラであれば速度の差はなくなる。
0873１３２人目の素数さん
垢版 |
2020/08/08(土) 23:06:27.95ID:UlZ+DmHh
一人だけエラー吐いて、まともに動かない関数書いた無能が、
それでも、恥ずかしげもなく、一生懸命マウント取ろうとあがく姿は、
滑稽で哀れになってくるわ
0874１３２人目の素数さん
垢版 |
2020/08/13(木) 17:34:36.29ID:5b54QocS
＞859がrec.prod(1:13)であっさりオーバーフローしてしまうのは、integer型だから
の意味が理解できなくて、実験してみた。

> rec.prod <- function(x) {
+ l <- length(x)
+ if (l == 0)
+ return(1)
+ x1.a <- abs(x[1])
+ if (l == 1)
+ x1.a
+ else
+ x1.a * rec.prod(x[-1])
+ }
>
>
> rec.prod(1:13)
[1] NA
Warning message:
In x1.a * rec.prod(x[-1]) : NAs produced by integer overflow
> rec.prod(as.numeric(1:13))
[1] 6227020800
> rec.prod(c(1,2,3,4,5,6,7,8,9,10,11,12,13))
[1] 6227020800
> rec.prod(c(1L,2L,3L,4L,5L,6L,7L,8L,9L,10L,11L,12L,13L))
[1] NA
Warning message:
In x1.a * rec.prod(x[-1]) : NAs produced by integer overflow

integer型のデータだとオーバーフローしちゃうけど、as.numericでdouble型にするとエラーがでない。

一方、

> f <- function(x){
+ if(length(x) == 0) return(1)
+ abs(x[1]) * f(x[-1])
+ }
> f(1:13)
[1] 6227020800
> f(as.numeric(1:13))
[1] 6227020800
> f(c(1,2,3,4,5,6,7,8,9,10,11,12,13))
[1] 6227020800
> f(c(1L,2L,3L,4L,5L,6L,7L,8L,9L,10L,11L,12L,13L))
[1] 6227020800
>

とエラーがでないなぁ。

何故だろう？
0875１３２人目の素数さん
垢版 |
2020/08/20(木) 13:26:56.75ID:mPUCNi08
Package ‘XXX’ was installed before R 4.0.0: please re-install it
とメッセージがでてXXXをインストールすると次は

Package ‘YYY’ was installed before R 4.0.0: please re-install it
とでてきてうんざりしていた。

stackoverflow.comで以下の答をみつけて解決した。

.libPaths()
update.packages(ask=FALSE, checkBuilt=TRUE)
0876１３２人目の素数さん
垢版 |
2020/09/09(水) 22:54:07.12ID:IR7822fG
ありがとうございます
0877１３２人目の素数さん
垢版 |
2020/10/19(月) 07:15:12.45ID:21zikI9w
PならばQ　をパイプで　P %=>%　Q　として論理演算で遊んでみた。

"
問題　：　「シリツ医 ならば (馬鹿 ならば 裏口 である)」という命題と同値な命題はどれか？

1 : シリツ医 ならば (裏口 ならば 馬鹿 である)
2 : 馬鹿 ならば (シリツ医 ならば 裏口 である)
3 : 馬鹿 ならば (裏口 ならば シリツ医 である)
4 : 裏口 ならば (シリツ医 ならば 馬鹿 である)
5 : 裏口 ならば (馬鹿 ならば シリツ医 である)
"

# PならばQ ≡　(P かつ (Qでない))ではない
'%=>%' = function(P,Q) !(P & !Q)
# premise
A = function(S,B,U) S %=>% (B %=>% U)

# choice
B1 = function(S,B,U) S %=>% (U %=>% B)
B2 = function(S,B,U) B %=>% (S %=>% U)
B3 = function(S,B,U) B %=>% (U %=>% S)
B4 = function(S,B,U) U %=>% (S %=>% B)
B5 = function(S,B,U) U %=>% (B %=>% S)

# premise => choice
C1 = function(S,B,U) A(S,B,U) %=>% B1(S,B,U)
C2 = function(S,B,U) A(S,B,U) %=>% B2(S,B,U)
C3 = function(S,B,U) A(S,B,U) %=>% B3(S,B,U)
C4 = function(S,B,U) A(S,B,U) %=>% B4 (S,B,U)
C5 = function(S,B,U) A(S,B,U) %=>% B5 (S,B,U)

# combination grid of TRUE & FALSE
gr=expand.grid(c(T,F),c(T,F),c(T,F))

# test for premise => choice
all(mapply(C1,gr[,1],gr[,2],gr[,3]))
all(mapply(C2,gr[,1],gr[,2],gr[,3]))
all(mapply(C3,gr[,1],gr[,2],gr[,3]))
all(mapply(C4,gr[,1],gr[,2],gr[,3]))
all(mapply(C5,gr[,1],gr[,2],gr[,3]))
■ このスレッドは過去ログ倉庫に格納されています
ニューススポーツなんでも実況